Vorwort

Liebe Studierende,

herzlich willkommen zur Vorlesung “Programmierkurs für Chemiker”. In diesem Kurs werden Sie die Grundlagen des Programmierens mit der Programmiersprache Python erlernen und viele nützliche numerische Algorithmen für die Chemie und darüber hinaus implementieren.

Diese Veranstaltung findet

• dienstags von 10:15 bis 11:45 Uhr und
• donnerstags von 14:15 bis 15:45 Uhr

im Seminarraum 00.006 des Instituts für Theoretische Chemie (Campus Hubland Nord) statt.

Syllabus

Ein vorläufiger Syllabus ist unten dargestellt. Dieser wird im Laufe des Semesters entsprechend dem Fortschritt der Vorlesung kontinuierlich aktualisiert.

Motivation

Warum sollte ich im Chemie-Studium programmieren lernen?

Im 21. Jahrhundert wird Programmieren zu einer immer wichtigeren Fähigkeit. Dies gilt nicht nur für die Informatik, sondern auch für die Chemie, Physik, Biologie und viele andere Wissenschaften. Die Kenntnis von Programmiersprachen kann helfen, eine Vielzahl von Problemen effizienter zu lösen und repetitive Aufgaben zu automatisieren. Darüber hinaus können Programmiersprachen verwendet werden, um Daten zu visualisieren und zu analysieren, was ein wichtiger Bestandteil wissenschaftlicher Arbeit ist.

Warum Python?

Python ist eine der am weitesten verbreiteten Programmiersprachen und besonders anfängerfreundlich. Die folgende Abbildung zeigt die Beliebtheit einiger Programmiersprachen gemäß dem PopularitY of Programming Language (PYPL) Index (Stand Jan. 2024).

Beliebtheit von Programmiersprachen. Entnommen aus dem PYPL Index.

Ein Grund dafür ist, dass Python oft deutlich weniger Code als kompilierte Sprachen wie Java oder C/C++ benötigt, um die gleichen Algorithmen zu implementieren.

Programmlänge, gemessen in der Anzahl der nicht-kommentierten Codezeilen (LOC).

Darüber hinaus ist es oft möglich, die gleiche Aufgabe mit einer Skriptsprache wie Python in deutlich kürzerer Zeit zu lösen als mit einer kompilierten Sprache.

Entwicklungszeit, um eine bestimmte Programmieraufgabe zu lösen, gemessen in Stunden.

Auf der anderen Seite argumentieren einige Leute, dass Python eine langsame Sprache ist. Dies ist bis zu einem gewissen Grad wahr. Im folgenden betrachten wir einen Grund, warum Python als langsam angesehen wird, aber auch, warum Python trotzdem sehr performant sein kann.

Python ist (auch) eine interpretierte Sprache

Python selbst ist ein C-Programm, das den Quellcode zuerst in sogenannten Bytecode kompiliert und diesen dann interpretiert und ausführt. Dies steht im Gegensatz zu kompilierten Sprachen wie C, C++, Rust, etc., bei denen der Quellcode in Maschinencode kompiliert wird. Der Compiler kann viele Optimierungen am Code vornehmen, was zu einer kürzeren Laufzeit führt.

s = 0
for i in range(100_000_000):
    if i % 2 == 1:
        s += i

#![allow(unused)]
fn main() {
let mut s: usize = 0;
for i in 0..100_000_000 {
    if i % 2 == 1 {
        s += i;
    }
}
}

Python kann sehr performant sein

Während dieser Veranstaltung werden wir oft Python-Bibliotheken wie NumPy oder SciPy für mathematische Algorithmen und insbesondere lineare Algebra verwenden. Diese Pakete bringen zwei wesentliche Vorteile. Einerseits ermöglichen sie die sehr einfache Verwendung komplizierter Algorithmen und andererseits sind diese Pakete in kompilierten Sprachen wie C oder Fortran geschrieben. Auf diese Weise können wir von den Leistungsvorteilen profitieren, ohne eine potenziell kompliziertere Sprache lernen zu müssen.

Erste Schritte

Python-Distribution

Wir empfehlen die Verwendung der Miniforge Distribution. Falls Sie bereits die Anaconda Distribution installiert haben, können Sie auch diese für den Kurs verwenden. Folgen Sie die untenstehenden Anweisungen, um Miniforge zu installieren. Sie sollten eine Python-Version von 3.9 oder höher verwenden.

Installation

1. Laden Sie den Miniforge-Installer gemäß Ihrem Betriebssystem hier herunter.
2. Führen Sie den Installer aus:
   • Mac OS & Linux: Öffnen Sie das Terminal, navigieren Sie zu Ihrem Downloads-Ordner und rufen Sie

     bash [Miniforge-Installer].sh
     

     auf, wobei [Miniforge-Installer] durch den Namen des Installers ersetzt werden sollte, z.B. Miniforge3-Linux-x86_64.sh oder Miniforge3-MacOSX-arm64.sh.
   • Windows: Doppelklicken Sie auf den Installer im Explorer.

Python-Pakete

Während der Vorlesung werden wir mehrere Python-Pakete verwenden, die bequem mit dem Paketmanager mamba installiert werden können, der mit der Miniforge-Distribution installiert wird. Einige der wichtigsten Pakete sind in der Tabelle unten aufgeführt. Bitte installieren Sie diese Pakete, indem Sie die folgenden Befehle in Ihrem Terminal ausführen.

Jupyter Notebook

Wenn Sie ein Anfänger in Python sind, empfehlen wir die Verwendung von Jupyter Notebook zum Schreiben und Ausführen Ihres Python-Codes. Es handelt sich um ein blockweise ausführbares Dokument, das Code, Text und Grafiken enthalten kann. Um Ihre Python-Installation zu testen und sich mit Jupyter Notebook vertraut zu machen, laden Sie bitte die Datei Crashkurs: Jupyter Notebooks und Python herunter, die in WueCampus bereitgestellt ist. Sie können das Notebook öffnen indem Sie

jupyter notebook

in Ihrem Terminal eingeben. Anschließend öffnet sich ein neuer Tab in Ihrem Webbrowser, wo Sie zum Ordner navigieren können, in dem sich das Jupyter Notebook befindet. Öffnen Sie das Notebook und führen Sie alle Code-Zellen aus. Wenn Ihre Installation erfolgreich war, sollten Sie keinen Fehler erhalten. Zur Bedienung des Notebooks finden Sie viele weitere Tutorials im Internet, wie z.B. hier.

Empfehlungen für weitere integrierte Entwicklungsumgebungen (IDEs) oder Editoren

Wenn Sie bereits mit Python vertraut sind und eine erweiterte Entwicklungsumgebung nutzen möchten, sollten Sie eine IDE oder einen spezialisierten Editor verwenden. Obwohl es möglich ist, Python-Code mit einem herkömmlichen Texteditor zu schreiben, kann eine gute IDE Ihre Programmiererfahrung erheblich verbessern. Deshalb stellen wir Ihnen hier einige Empfehlungen für IDEs und spezialisierte Editoren neben Jupyter Notebook vor.

• Spyder

  • voll ausgestattete Python-IDE mit Schwerpunkt auf wissenschaftlichr Entwicklung
  • leichtgewichtig und einfach zu konfigurieren

• Visual Studio Code

  • einer der meistgenutzten Editoren trotz seiner Leichtigkeit
  • fast alle Funktionen einer IDE, obwohl es offiziell keine ist

• PyCharm

  • kommerzielle (für Studierende kostenlos) IDE mit vielen Funktionen
  • für sehr große und komplexe Python-Projekte geeignet
  • möglicherweise zu umständlich für kleine Projekte wie die in dieser Veranstaltung

• Vim/NeoVim

  • Kommandozeilen-Editor, der auf fast allen Unix-ähnlichen Computern vorinstalliert ist
  • möglicherweise sehr unhandlich am Anfang aufgrund der vielen Tastenkombinationen und der mangelnden Anfängerfreundlichkeit im Vergleich zu einer typischen IDE oder Jupyter Notebooks
  • extrem konfigurierbar
  • mit einigem Aufwand zu einem sehr umfangreichen und komfortablen Editor anpassbar
  • immer noch einer der meistgenutzten Editoren in der Softwareentwicklung

Bedienung dieser Webseite

Dieser Abschnitt gibt eine Einführung in die Bedienung dieses Vorlesungsskripts. Dieses ist in Kapitel organisiert, welche wiederum in Abschnitte unterteilt sind. Typischerweise ist jedes (Unter-)Kapitel in eine Reihe von Überschriften unterteilt.

Navigation

Es gibt mehrere Methoden, um durch die Kapitel eines Buches zu navigieren.

Die Seitenleiste auf der linken Seite bietet eine Liste aller Kapitel. Durch Klicken auf einen der Kapiteltitel wird die entsprechende Seite geladen.

Die Seitenleiste erscheint möglicherweise nicht automatisch, wenn das Fenster zu schmal ist, insbesondere auf mobilen Displays. In diesem Fall kann das Menüsymbol (drei horizontale Balken) oben links auf der Seite gedrückt werden, um die Seitenleiste zu öffnen und zu schließen.

Die Pfeiltasten links und rechts neben der Seite können verwendet werden, um zum vorherigen oder nächsten Kapitel zu navigieren.

Die Pfeiltasten auf der Tastatur können verwendet werden, um zum vorherigen oder nächsten Kapitel zu navigieren.

Obere Menüleiste

Die Menüleiste oben auf der Seite bietet einige Symbole zur Interaktion mit dem Vorlesungsskript.

Mit einem Klick auf die Menüleiste scrollt die Seite nach oben.

Suche

Dieses Vorlesungsskript verfügt über ein eingebautes Suchsystem. Durch Drücken des Suchsymbols () in der Menüleiste oder Drücken der Taste S auf der Tastatur wird ein Eingabefeld zum Eingeben von Suchbegriffen geöffnet. Das Eingeben von Begriffen zeigt passende Kapitel und Abschnitte in Echtzeit an.

Durch Klicken auf eines der Ergebnisse wird zu diesem Abschnitt gesprungen. Die Pfeil-nach-oben- und Pfeil-nach-unten-Tasten können verwendet werden, um die Ergebnisse zu navigieren, und Enter öffnet den markierten Abschnitt.

Nach dem Laden eines Suchergebnisses werden die passenden Suchbegriffe im Text hervorgehoben. Durch Klicken auf ein hervorgehobenes Wort oder Drücken der Esc-Taste werden die Hervorhebungen entfernt.

Codeblöcke

Codeblöcke besitzen ein Kopieren-Symbol , das den Codeblock in die lokale Zwischenablage kopiert.

Hier ein Beispiel:

print("Hello, World!")

Wir verwenden häufig die assert-Anweisung in Codeblöcken, um Ihnen den Wert einer Variablen anzuzeigen. Da die Codeblöcke in diesem Dokument nicht interaktiv sind (Sie können sie nicht einfach in Ihrem Browser ausführen), ist es nicht möglich, den Wert der Variablen auf dem Bildschirm auszugeben. Deshalb stellen wir für Sie sicher, dass alle Codeblöcke in diesem Vorlesungsskript fehlerfrei laufen. Der folgende Codeblock zeigt Ihnen beispielsweise an, dass der Variable „a“ der Wert 2 zugewiesen wurde:

a = 2
assert a == 2

Regressionsanalyse

Eine sehr häufige Aufgabe der Datenanalyse besteht darin, ein Modell an einen Datensatz anzupassen. Dies nennt man Regressionsanalyse, ein grundlegendes statistisches Werkzeug in unzähligen wissenschaftlichen Untersuchungen, auch in der Chemie, welches die Beziehungen zwischen Variablen aufdeckt.

Betrachten Sie die folgenden praktischen Beispiele: Anwendung des Lambert- Beerschen Gesetzes zur Quantifizierung der Konzentration einer Lösung basierend auf Lichtabsorption, Bestimmung des Einflusses der Temperatur auf Reaktionsgeschwindigkeiten, Herausfinden der Beziehung zwischen Substratkonzentration und Enzymaktivität. Jedes dieser Beispiele zeigt, wie die Regressionsanalyse dabei hilft, experimentelle Daten zu modellieren, Beziehungen zwischen Variablen zu extrahieren und Vorhersagen zu treffen.

In diesem Kapitel werden wir uns mit den Grundlagen der Regressionsanalyse befassen, angefangen bei der einfachen linearen Regression bis hin zu komplexeren Modellen. Wir werden uns auch mit den zugrundeliegenden Optimierungsalgorithmen befassen, die zum Fitten dieser Modelle an die Daten verwendet werden. Darüberhinaus werden wir die Bedeutung der Modellevaluation und der Regularisierungstechniken zur Vermeidung vom Überfitten diskutieren.

Methode der kleinsten Quadrate

Das Ziel aller Regressionsprobleme besteht darin, ein Modell zu finden, welches am besten zu den Daten passt und hoffentlich zur Vorhersage neuer Datenpunkte verwendet werden kann. Einfachheitshalber gehen wir davon aus, dass unsere Daten nur eine unabhängige Variable ($x$) und eine abhängige Variable ($y$) haben. Eine Verallgemeinerung auf mehrere unabhängige Variablen ist simpel und wird in einem späteren Kapitel diskutiert. Darüber hinaus gehen wir davon aus, dass die Daten reellwertig sind. Für insgesamt $N$ Datenpunkte schreiben wir $$\begin{array}{rl}x& =(x_1,x_2,\dots ,x_N)\in {\mathbb{R}}^N\\ y& =(y_1,y_2,\dots ,y_N)\in {\mathbb{R}}^N,\end{array}$$ d.h. $x$ und $y$ sind $N$-dimensionale Vektoren mit reellwertigen Komponenten $x_i$ und $y_i$.

Wir können die Beziehung zwischen $x_i$ und $y_i$ als eine Funktion $f$ betrachten. Unter idealisierten Bedingungen haben wir $$y_i=f(x_i)+{ϵ}_i,\text{\hspace{1em}(1.1)}$$ wobei die Funktion $f$ die wahre Beziehung zwischen den Datenpunkten $x_i$ und $y_i$ beschreibt, während ${ϵ}_i$ den zufälligen statistischen Fehler darstellt.

In der realen Welt ist die Beziehung $f$ jedoch oft unbekannt. Daher müssen wir ein Modell verwenden, um sie abzuschätzen. In diesem Fall schreiben wir $$y_i=\hat{f}(\beta ;x_i)+{\widehat{ϵ}}_i,\text{\hspace{1em}(1.2)}$$ wobei $\hat{f}(\beta ;x_i)$ der Abschätzer des Modells von $f$ ist, parametrisiert durch $\beta \in {\mathbb{R}}^n$ mit $n$ Parametern im Modell. Der Fehlerterm ${\widehat{ϵ}}_i$ ist nun eine Kombination aus dem statistischen Fehler und dem unmodellierten Teil der echten Beziehung. Das Ziel der Regressionsanalyse besteht darin, die Parameter $\beta$ zu finden, sodass das Modell die Daten am besten wiedergibt.

Aber was bedeutet “am besten”? Ein üblicher Ansatz ist die Methode der kleinsten Quadrate, die darauf abzielt, die Summe der quadratischen Fehler ${\widehat{ϵ}}_i$ zu minimieren, d.h. die Verlustfunktion der kleinsten Quadrate $$L(\beta ;x,y)=\sum _{i=1}^N{\widehat{ϵ}}_i^2=\sum _{i=1}^N(y_i-\hat{f}(\beta ;x_i){)}^2.\text{\hspace{1em}(1.3)}$$ zu minimieren. Mathematisch ausgedrückt wollen wir die Parameter ${\beta }^{\ast }$ durch Lösen des Optimierungsproblems $${\beta }^{\ast }=\underset{\beta \in {\mathbb{R}}^n}{\mathrm{argmin}}L(\beta ;x,y)=\underset{\beta \in {\mathbb{R}}^n}{\mathrm{argmin}}\sum _{i=1}^N(y_i-\hat{f}(\beta ;x_i){)}^2.\text{\hspace{1em}(1.4)}$$ finden. Wenn wir den Vektor der vorhergesagten Werte $\hat{y}$ als $${\hat{y}}_i=\hat{f}(\beta ;x_i),$$ schreiben, können wir Gl. (1.4) als $${\beta }^{\ast }=\underset{\beta \in {\mathbb{R}}^n}{\mathrm{argmin}}\sum _{i=1}^N(y_i-{\hat{y}}_i{)}^2=\underset{\beta \in {\mathbb{R}}^n}{\mathrm{argmin}}\Vert y-\hat{y}{\Vert }_2^2,\text{\hspace{1em}(1.5)}$$ formulieren, wobei $\Vert v{\Vert }_2$ die euklidische Norm oder die ${\ell }_2$-Norm eines Vektors $v$ bezeichnet, definiert als $$\Vert v{\Vert }_2=\sqrt{\sum _{i=1}^N{v}_i^2}.$$

Die Methode der kleinsten Quadrate ist eine beliebte Wahl für die Regressionsanalyse, da sie geschlossene Lösungen für einige einfache, aber wichtige Modelle bietet und mit Methoden der numerischen linearen Algebra für komplexere Modelle effizient gelöst werden kann. Es weist jedoch einige Nachteile auf, z.B. die Empfindlichkeit gegenüber Ausreißern und die Anfälligkeit für Überanpassungen. Eine Alternative zur Methode der kleinsten Quadrate ist die Methode der kleinsten absoluten Abweichungen, welche die ${\ell }_1$-Norm anstelle der ${\ell }_2$-Norm verwendet. Diese ist für einen Vektor $v$ definiert als $$\Vert v{\Vert }_1=\sum _{i=1}^N|v_i|.$$

Das Optimierungsproblem wird dann formuliert als $${\beta }^{\ast }=\underset{\beta \in {\mathbb{R}}^n}{\mathrm{argmin}}\Vert y-\hat{y}{\Vert }_1.\text{\hspace{1em}(1.6)}$$

Es gibt viele andere Verlustfunktionen, die für die Regression verwendet werden können, z.B. die Huber-Verlustfunktion, die log-cosh-Verlustfunktion, die Quantilverlustfunktion usw. Obwohl die Methode der kleinsten Quadrate in vielen Fällen gut funktioniert, ist es wichtig, sich der Grenzen der Methode bewusst zu sein und gegebenenfalls andere Verlustfunktionen heranzuziehen.

Lineare Regression

Eines der einfachsten Modell der Regressionsanalyse ist das lineare Modell, wobei der Abschätzer $\hat{f}(\beta ;x_i)$ durch eine lineare Funktion der Form $$\hat{f}(\beta ;x_i)={\beta }_0+{\beta }_1{x}_i\text{\hspace{1em}(1.7)}$$ mit den skalaren Parametern ${\beta }_0$ und ${\beta }_1$ gegeben ist. Die Regressionsanalyse mit dem linearen Modell wird als lineare Regression bezeichnet.

Theoretische Grundlagen

Setzen wir nun das lineare Modell nach Gl. (1.7) in die Verlustfunktion der kleinsten Quadrate gemäß Gl. (1.3) ein, erhalten wir $$L(\beta ;x,y)=\sum _{i=1}^N(y_i-{\beta }_0-{\beta }_1{x}_i{)}^2.\text{\hspace{1em}(1.8)}$$

Da unser Ziel darin besteht, die Parameter $\beta$ zu finden, welche die Verlustfunktion minimieren (vgl. Gl. (1.4)), müssen zumindest die notwendigen Bedingungen dafür erfüllt werden, d.h. die partiellen Ableitungen von $L$ nach ${\beta }_0$ und ${\beta }_1$ müssen verschwinden:

$$\frac{\partial L}{\partial {\beta }_0}\stackrel{!}{=}0\text{und}\frac{\partial L}{\partial {\beta }_1}\stackrel{!}{=}0.$$

Nach einer etwas länglichen aber einfachen

erhalten wir die Matrixgleichung für die Parameter ${\beta }_0$ und ${\beta }_1$: $$\underset{\mathbf{A}}{\underbrace{\left(\begin{array}{cc}N& \sum _{i=1}^N{x}_i\\ \sum _{i=1}^N{x}_i& \sum _{i=1}^N{x}_i^2\end{array}\right)}}\underset{\vec{x}}{\underbrace{\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right)}}=\underset{\vec{b}}{\underbrace{\left(\begin{array}{c}\sum _{i=1}^N{y}_i\\ \sum _{i=1}^N{x}_i{y}_i\end{array}\right)}}.\text{\hspace{1em}(1.9)}$$ Die Lösung dieses Gleichungssystems liefert uns die optimalen Parameter ${\beta }_0$ und ${\beta }_1$ für das lineare Modell.

Die (zumindest formale) Lösung des Gleichungssystems (1.9) ist $\vec{x}={\mathbf{A}}^{-1}\vec{b}$, mit der Inversen der Matrix $\mathbf{A}$. Eine Matrix ist genau dann invertierbar, wenn ihre Determinante ungleich null ist. Die Determinante der Systemmatrix $\mathbf{A}$ ist $$\begin{array}{rl}\det (\mathbf{A})& =N\sum _{i=1}^N{x}_i^2-{\left(\sum _{i=1}^N{x}_i\right)}^2\\ & =N^2\left[\frac{1}{N}\sum _{i=1}^N{x}_i^2-{\left(\frac{1}{N}\sum _{i=1}^N{x}_i\right)}^2\right]\\ & =N^2\left[\frac{1}{N}\sum _{i=1}^N{x}_i^2-{\bar{x}}^2\right]\\ & =N^2\left[\frac{1}{N}\sum _{i=1}^N(x_i-\bar{x}{)}^2\right],\end{array}$$ wobei wir die Abkürzung $\bar{x}=\frac{1}{N}{\sum }_{i=1}^N{x}_i$ für den Mittelwert der $x_i$ eingeführt haben. Hätten wir nur einen einzigen Datenpunkt, wäre der Mittelwert gleich dem einzigen Datenpunkt und die Determinante der Systemmatrix gleich null. Das bedeutet, dass die Matrix $\mathbf{A}$ nicht invertierbar ist und wir keine eindeutige Lösung für das lineare Modell erhalten. Erst ab zwei verschiedenen Datenpunkten ist das Gleichungssystem (1.9) eindeutig lösbar.

Nun wollen wir mit Hilfe der Gl. (1.9) die lineare Regression für ein einfaches Beispiel implementieren.

Implementierung

Betrachten wir die folgenden Messdaten¹ der Lambert-Beer-Beziehung für Methylenblau in Wasser bei verschiedenen Konzentrationen $c$ und den zugehörigen Absorbanzen $A$, gemessen bei $610\mathrm{nm}$ mit einer Schichtdicke $d$ von 1 cm:

Bevor wir fortfahren können, müssen wir die Daten irgendwie in Python importieren. Der einfachste Weg für so einen kleinen Datensatz ist die manuelle Eingabe.

concentrations = [
    2.125, 4.250, 6.375, 8.500, 10.63, 12.75, 14.88, 17.00, 19.13, 21.25,
    23.38, 25.50, 27.63, 29.75, 31.88, 34.00, 36.13, 38.25, 40.38, 42.50,
]
absorbances = [
    0.0572, 0.1391, 0.2049, 0.2754, 0.3420, 
    0.4139, 0.4956, 0.5815, 0.6806, 0.7481,
    0.8242, 0.9130, 1.0043, 1.0809, 1.1511,
    1.2483, 1.3373, 1.4027, 1.4927, 1.5853,
]

Hier haben wir die Daten in den Variablen concentrations und absorbances der Typ List (Liste) definiert. Das erkennt man an der Verwendung der eckigen Klammern [] und dem Komma , zwischen den einzelnen Werten. Mit dem Gleichheitszeichen = weisen wir den Variablen die Werte zu.

Obwohl Python standardmäßig schon einige mathematische Funktionen für Listen bereitstellt, hat das Paket numpy noch mehr sehr nützliche Operationen für solchen Datenstrukturen. Daher importieren wir numpy und konvertieren die Listen in den Datentyp numpy-array.

import numpy as np

concentrations = np.array(concentrations)
absorbances = np.array(absorbances)

Da numpy etwas länger zu schreiben ist, verwenden wir den Alias np, definiert durch import numpy as np. Ab sofort können wir Inhalte dieses Pakets mit np ansprechen. In den folgenden Zeilen verwenden wir die Funktion np.array um die Listen in Arrays zu konvertieren. Hier wird es deutlich, dass das Gleichheitszeichen = beim Programmieren eine andere Bedeutung hat als in der Mathematik. Während in der Mathematik das Gleichheitszeichen eine Äquivalenz zwischen zwei Seiten ausdrückt, stellt dieses in der Programmierung eine Zuweisung dar. Deshalb ist es hier möglich, die gleiche Variablename für die Arrays zu verwenden, da die rechte Seite des Gleichheitszeichens zuerst ausgewertet wird und dann der Wert der rechten Seite der Variablen auf der linken Seite zugewiesen wird.

Als nächstes wollen wir die Elemente von $\mathbf{A}$ und $\vec{b}$ in Gl. (1.9) berechnen. Dafür stellen wir erstmal sicher, dass unsere Datenarrays gleich lang sind. Danach können wir die Elemente mit Hilfe der np.sum-Funktion berechnen.

# make sure the number of data points is the same
assert len(concentrations) == len(absorbances)

number_of_points = len(concentrations)
sum_x = np.sum(concentrations)
sum_x_sq = np.sum(concentrations**2)
sum_y = np.sum(absorbances)
sum_xy = np.sum(concentrations * absorbances)

Hier haben wir die built-in Funktion len verwendet, um die Länge der Arrays zu bestimmen. Die Funktion np.sum berechnet die Summe aller Elemente in einem Array. Mit dem Stern * ist die Multiplikation gemeint. Die Verwendung von * zwischen zwei Arrays führt zu einer elementweisen Multiplikation. Der Doppelstern ** bedeutet in Python die Potenzierung. Hier wird also das Array concentrations elementweise quadriert.

Nun können wir die Systemmatrix $\mathbf{A}$ und den Vektor $\vec{b}$ zusammenstellen:

a_arr = np.array([
    [number_of_points, sum_x],
    [sum_x, sum_x_sq],
])
b_arr = np.array([sum_y, sum_xy])

Beachte Sie, dass für die Erstellung einer Matrix (bzw. eines sog. 2D-Arrays) eine Liste von Listen verwendet wird. Die innere Liste entspricht einer Zeile der Matrix. Die äußere Liste enthält die Zeilen der Matrix. Jetzt können wir endlich das Gleichungssystem (1.9) lösen, um die optimalen Parameter ${\beta }_0$ und ${\beta }_1$ zu erhalten:

beta = np.linalg.solve(a_arr, b_arr)
print(beta)

Die Funktion np.linalg.solve löst (numerisch) ein lineares Gleichungssystem. Mit der print-Funktion geben wir die Lösung anschließend aus.

Wenn man schon weiß, welche Werte ${\beta }_0$ und ${\beta }_1$ haben sollen und nur den Algorithmus testen möchte, kann man die Verifikation z.B. wie folgt durchführen:

beta0 = beta[0]
beta1 = beta[1]
assert np.isclose(beta0, -0.04907034)
assert np.isclose(beta1, 0.03800109)

Hier haben wir zuerst den Array beta mit eckigen Klammern [] indexiert. Beachten Sie, dass die Indizierung in Python bei 0 beginnt. Das bedeutet, dass ${\beta }_0$ dem 0-ten Eintrag und ${\beta }_1$ den ersten Eintrag des Arrays beta entspricht. Anschließend haben wir die Werte von ${\beta }_0$ und ${\beta }_1$ mit Referenzwerten verglichen.

Das Lambert-Beer-Gesetz besagt, dass die Absorption $A$ linear von der Konzentration $c$ abhängt mit der Proportionalitätskonstante $\epsilon d$, also $$A=\epsilon dc.$$ Das molare Extinktionskoeffizient $\epsilon$ kann also aus dem linearen Parameter ${\beta }_1$ berechnet werden als $$\epsilon ={\beta }_1/d=0.038{\mu \mathrm{M}}^{-1}\cdot \mathrm{c}{\mathrm{m}}^{-1}=38000{\mathrm{M}}^{-1}\cdot \mathrm{c}{\mathrm{m}}^{-1}.$$

Jetzt haben wir die optimalen Parameter ${\beta }_0$ und ${\beta }_1$ für den obigen Datensatz berechnet. Wie soll man wissen, wie gut die lineare Regression ist? Dafür gibt es verschiedene mathematische Gütekriterien, die wir hier aber erstmal nicht betrachten. Stattdessen wollen wir die optimalen Parameter grafisch darstellen und die Qualität der Regression visuell beurteilen.

Visualisierung

Zur Visualisierung der optimalen linearen Regression verwenden wir das Python-Paket matplotlib. Hier importieren wir das Untermodul pyplot mit dem Alias plt:

import matplotlib.pyplot as plt

Für die graphische Darstellung mit matplotlib benötigen wir immer die Objekte Figure und Axes. Diese können wir mit der Funktion plt.subplots erstellen:

fig, ax = plt.subplots(figsize=(8, 6))

Wir haben hier als Argument figsize=(8, 6) übergeben, welches die Größe der Abbildung in Zoll angibt. Die Größe der Abbildung ist aber nicht in Stein gemeißelt und kann je nach Bedarf angepasst werden.

Nun wollen wir die Messdaten als Punkte und die lineare Regression als Linie in das Diagramm eintragen. Dafür verwenden wir die Funktion ax.plot:

ax.plot(concentrations, absorbances, 'o', label='data')
ax.plot(concentrations, beta0 + beta1 * concentrations, label='fit')

Für die Messdaten haben wir das Argument o verwendet, damit diese als Punkte dargestellt werden. Andere Marker wie z.B. s für Quadrate oder x für Kreuze können auch verwendet werden. Ohne Angabe dieses Arguments werden die Daten als Linie dargestellt. Bei beiden Aufrufen von ax.plot haben wir die Argumente label verwendet. Diese dienen dazu, die Linien in der Legende zu benennen.

Für die Vollständigkeit wollen wir auch die Achsen noch beschriften und eine Legende hinzufügen:

ax.set_xlabel('concentration / µM')
ax.set_ylabel('absorbance')

# automatically create a legend
ax.legend()

Die Funktionen ax.set_xlabel und ax.set_ylabel akzeptieren einen str (String, Zeichenkette) als Argument, der als Achsenbeschriftung dient. Die Funktion ax.legend fügt eine Legende hinzu. Ruft man sie ohne Argumente auf, wird die Legende aus den label-Argumenten der ax.plot-Funktionen automatisch generiert.

Nachdem alle Elemente des Diagramms eingefügt wurden, können wir die Abbildung mit der Funktion plt.show anzeigen lassen:

plt.show()

Wenn das Programm erfolgreich ausgeführt wird, sollte ein Diagramm wie das Folgende erscheinen:

Dieses Diagramm zeigt uns die Ergebnisse der linearen Regression, auch wenn man über die ästhetische Gestaltung des Diagramms streiten kann. Im Laufe dieses Kurses werden wir noch mehr Funktionalitäten von matplotlib kennenlernen, die uns helfen, Elemente des Diagramms dem persönlichen Geschmack entsprechend zu gestalten.

Auf den ersten Blick scheint die Regressionsgerade sehr gut zu den Messdaten zu passen. Um den Unterschied zwischen den Messdaten und der Regressionsgerade deutlicher zu machen, können wir die Residuen, also die Differenz zwischen den Messdaten und der Regressionsgerade, in einem weiteren Diagramm darstellen.

residuals = absorbances - (beta0 + beta1 * concentrations)

fig2, ax2 = plt.subplots(figsize=(8, 6))
ax2.bar(concentrations, residuals)
ax2.set_xlabel('concentration / µM')
ax2.set_ylabel('absorbance residuals')
plt.show()

Nach dem Berechnen und Speichern der Residuen in der Variablen residuals wurden ein weiteres Figure- sowie Axes-Objekt erstellt. Die Residuen wurden in ein Balkendiagramm dargestellt, welches mit der Methode ax2.bar erstellt werden kann. Das Diagramm sieht dann wie folgt aus:

Jetzt kann man erkennen, dass bei niedrigen und hohen Konzentrationen die Abweichung positiv ist, während sie bei mittleren Konzentrationen negativ ist. Das könnte auf eine leichte positive Krümmung der Daten hinweisen, die mit dem linearen Modell nicht erfasst wird. Eine Zuordnung dieses Verhaltens zu einer systematischen oder zufälligen Abweichung bedarf allerdings in der Regel einer tieferen Analyse. Tatächlich liegt hier in den Daten eine leichte Nichtlinearität vor, welches die Abweichung vom Lambert-Beer-Gesetz v.a. bei hohen Konzentrationen demonstriert. Wer sich für die genaue Erklärung dieser Abweichung von Methylenblau in Wasser interessiert, kann sich z.B. die Publikation A. Fernández-Pérez, T. Valdés-Solís, G. Marbán, Dyes and Pigments 2019, 161, 448–456 durchlesen. Ob die lineare Regression in diesem Fall sinnvoll ist, hängt von der gewünschten Genauigkeit der Modellierung ab.

import numpy as np
import matplotlib.pyplot as plt

# Rest des Skripts

Übung

Aufgabe 1.1: Lineare und quadratische Regression

In dem obigen Beispiel haben wir die numpy Funktion np.linalg.solve verwendet, um die Lösung des Gleichungssystems der linearen Regression numerisch zu berechnen. In diesem Zusammenhang bedeutet dies, dass der Computer einer Reihe von Rechenschritten und Algorithmen folgt, um die approximative Lösung des Gleichungssystems zu finden. Für das Gleichungssystem der linearen Regression gibt es jedoch auch eine analytische Lösung, die direkt berechnet werden kann.

(a) Analytische Lösung der linearen Regression herleiten

Zeigen Sie, dass die Lösung des Gleichungssystems in Matrixform $$\underset{\mathbf{A}}{\underbrace{\left(\begin{array}{cc}N& \sum _{i=1}^N{x}_i\\ \sum _{i=1}^N{x}_i& \sum _{i=1}^N{x}_i^2\end{array}\right)}}\underset{\vec{x}}{\underbrace{\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right)}}=\underset{\vec{b}}{\underbrace{\left(\begin{array}{c}\sum _{i=1}^N{y}_i\\ \sum _{i=1}^N{x}_i{y}_i\end{array}\right)}}.$$ gegeben ist durch:

$$\begin{array}{rl}{\beta }_0& =\bar{y}-{\beta }_1\bar{x}\\ {\beta }_1& =\frac{{\sum }_{i=1}^N(x_i-\bar{x})(y_i-\bar{y})}{{\sum }_{i=1}^N(x_i-\bar{x}{)}^2}\end{array}$$

Lösen Sie dazu zunächst die erste Gleichung des Systems nach ${\beta }_0$ auf und setzen Sie das Ergebnis in die zweite Gleichung ein. Verwenden Sie außerdem die Definitionen der Mittelwerte $\bar{x}$ und $\bar{y}$.

(b) Implementieren der analytischen Lösung für Messdaten von Methylenblau

Nutzen Sie die analytische Lösung, um die Parameter der linearen Regression für die Messdaten von Methylenblau explizit zu berechnen. Vergleichen Sie die Ergebnisse mit den Ergebnissen, die Sie mit np.linalg.solve erhalten haben.

(c) Matrixgleichung der quadratischen Regression herleiten

Die quadratische Regression ist eine Erweiterung der linearen Regression, bei der die abhängige Variable $y$ durch ein Polynom zweiten Grades in der unabhängigen Variable $x$ angenähert wird. Die allgemeine Form der quadratischen Regression ist gegeben durch:

$$\hat{f}(\beta ;x_i)={\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2$$

In Analogie zur linearen Regression können wir die quadratische Regression als ein lineares Modell in den Parametern $\beta =({\beta }_0,{\beta }_1,{\beta }_2)$ auffassen. Zeigen Sie, dass dieses Modell durch die folgende Matrixgleichung beschrieben wird:

$$\left(\begin{array}{ccc}N& \sum _{i=1}^N{x}_i& \sum _{i=1}^N{x}_i^2\\ \sum _{i=1}^N{x}_i& \sum _{i=1}^N{x}_i^2& \sum _{i=1}^N{x}_i^3\\ \sum _{i=1}^N{x}_i^2& \sum _{i=1}^N{x}_i^3& \sum _{i=1}^N{x}_i^4\end{array}\right)\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\\ {\beta }_2\end{array}\right)=\left(\begin{array}{c}\sum _{i=1}^N{y}_i\\ \sum _{i=1}^N{x}_i{y}_i\\ \sum _{i=1}^N{x}_i^2{y}_i\end{array}\right)$$

Setzen Sie dazu die quadratische Funktion $\hat{f}(\beta ;x_i)$ in die allgemeine Form der Verlustfunktion der Methode der kleinsten Quadrate ein und bilden Sie die Ableitungen nach den gesuchten Parametern.

(d) Quadratische Regression implementieren und an Methylenblau-Daten anwenden

Fahren Sie nun fort wie für die lineare Regression, indem Sie das Gleichungssystem der quadratischen Regression aus Teilaufgabe (c) für die Methylenblau-Daten numerisch lösen. Konstruieren Sie dazu zunächst die benötigte Matrix, bzw. den Vektor in Form von Arrays, und verwenden Sie die Funktion np.linalg.solve. Plotten Sie anschließend die quadratische Regression zusammen mit den Datenpunkten. Plotten Sie ebenfalls die Resiuduen und vergleichen Sie die Ergebnisse mit der linearen Regression.

Numerische Optimierung

Die numerische Optimierung bietet uns die Möglichkeit, komplexe Funktionen, wie die Verlustfunktionen der kleinsten Quadrate (Gl. (1.3)), aber auch Funktionen im anderen Kontext, wie z.B. die Energie eines Moleküls, zu minimieren oder maximieren. Da die Maximierung einer Funktion $f$ äquivalent zur Minimierung von $-f$ ist, werden wir im Folgenden nur noch vom Minimieren sprechen. Die Fähigkeit, Optimierungsprobleme anzugehen, für die keine geschlossenen Lösungsformeln existieren, oder die Auswertung des analytischen Ausdrucks zu aufwendig ist, erweitert unser Werkzeugset in der Datenanalyse signifikant. Insbesondere erlaubt sie es, Lösungen für Modelle zu finden, die durch Nichtlinearitäten, hohe Dimensionalitäten oder ungewöhnliche Datenverteilungen charakterisiert sind.

Theoretische Grundlagen

Ein besonders zugänglicher und grundlegender Ansatz für die numerische Optimierung ist das Gradientenverfahren (engl. Gradient Descent). Dies ist ein iteratives Verfahren, welches von einem gegebenen Startpunkt $x^0$ ausgeht und in jedem Schritt der Richtung des steilsten Abstiegs der Funktion folgt. Mathematisch formuliert heißt das $$x^{k+1}=x^k-\alpha \nabla f(x^k)\text{\hspace{1em}(1.10)}$$ wobei $x^k$ der Schätzwert des Minimums im $k$-ten Schritt ist. Das hochgestellte $k$ hat hier nichts mit Potenzierung zu tun, sondern ist lediglich eine Notation, um die Verwechselung mit der $i$-ten Komponente des Vektors, die hier tiefgestellt wird, zu vermeiden. Der Gradient der Objektivfunktion $f$ bei $x^k$ kann dann als $\nabla f(x^k)={\left(\frac{\partial f}{\partial x_1^k},\dots ,\frac{\partial f}{\partial x_n^k}\right)}^{⊺}$ notiert werden. Die Propotionialitätskonstante $\alpha$ wird als Schrittweite oder learning rate bezeichnet. Das Verfahren wird solange wiederholt, bis eine oder mehrere Abbruchbedingungen erfüllt sind. Typische Abbruchbedingungen für iterative Optimierungsverfahren sind:

• Die Änderung des Funktionswertes ist kleiner als ein Schwellenwert
• Die Änderung des Schätzwertes ist kleiner als ein Schwellenwert
• Eine maximale Anzahl an Iterationen ist erreicht
• Die Norm des Gradienten ist kleiner als ein Schwellenwert.

Die Schrittweite $\alpha$ ist der einzige und zugleich ein wichtiger Parameter des Gradientenverfahrens, da sie die Konvergenzgeschwindigkeit und die Stabilität des Verfahrens beeinflusst. Ein zu kleiner Wert für $\alpha$ kann dazu führen, dass das Verfahren sehr langsam konvergiert, während ein zu großer Wert dazu führen kann, dass das Verfahren divergiert.

Damit wir das Gradientenverfahren nutzen können, müssen wir Zugang zum Gradienten der Objektivfunktion haben. Liegt dieser nicht analytisch vor, muss eine numerische Approximation verwendet werden. Ein einfacher Ansatz liefert die Methode der Finite Differenz . Hierbei wird die tangente der partiellen Ableitung durch die Sekante ersetzt, was zu einer Approximation der Form $$\frac{\partial f}{\partial x_i^k}\approx \frac{f(x^k+h{\hat{e}}_i)-f(x^k-h{\hat{e}}_i)}{2h}\text{\hspace{1em}(1.11)}$$ führt. Zudem ist ${\hat{e}}_i$ der $i$-te Einheitsvektor und $h$ ein kleiner Wert, der die Schrittweite der Approximation bestimmt. Genau genommen stellt Gl. (1.11) die zentrale finite Differenz 2. Ordnung dar. Es gibt auch einseitige Approximationen und Approximationen höherer Ordnung, die wir hier aber nicht weiter betrachten.

Implementierung

Finite Differenz

Als erstes implementieren wir die finite Differenz. Da wir mehrmals Ableitungen berechnen werden müssen, ist es sinnvoll, eine Funktion zu implementieren. Funktionen im Programmierkontext sind ähnlich zu mathematischen Funktionen, die eine Eingabe in eine Ausgabe umwandeln. Sie können aber noch einiges mehr.

Wir importieren zuerst numpy und die Objekte Callable und Any aus dem Modul typing:

import numpy as np
from typing import Callable, Any

Dann definieren wir die Funktion finite_difference, die den Gradienten einer Funktion func an der Stelle x0 berechnet.

def finite_difference(
    func: Callable[[np.ndarray, Any], float],
    x0: np.ndarray,
    h: float = 1e-5,
    args: tuple = (),
) -> np.ndarray:
    n: int = len(x0)
    grad: np.ndarray = np.zeros(n)

    for i in range(0, n):
        e: np.ndarray = np.zeros(n)
        e[i] = 1
        grad[i] = (
            func(x0 + h * e, *args) - func(x0 - h * e, *args)
        ) / (2 * h)

    return grad

Der erste Block der Funktion ist die Signatur, die angibt, welche Argumente die Funktion erwartet und welchen Typ die Rückgabe hat. Hierfür wurden die Objekte Callable und Any benutzt. Obwohl die explizite Angabe der Datentypen in Python optional und immer noch selten ist, ist es nach Meinung der Autoren eine gute Praxis. Näheres dazu finden Sie in dem folgenden

#![allow(unused)]
fn main() {
fn powi(x: f64, n: i32) -> f64 {
    let mut result: f64 = 1.0;
    for _ in 0..n {
        result *= x;
    }
    result
}
}

def powi(x, n):
    result = 1.0
    for _ in range(n):
        result *= x
    return result

Die Signatur der Funktion finite_difference besagt, dass sie eine Funktion func erwartet, die als Argument einen np.ndarray und irgendwas (Any) akzeptiert und einen Float zurückgibt. Die weiteren Argumente sind der Punkt x0 vom Typ np.ndarray, an dem der Gradient berechnet werden soll, die Schrittweite h vom Typ Float und ein Tuple (engl. tuple) args, das zusätzliche Argumente für die Funktion func liefern kann. Außerdem haben wir in der Signatur voreingestellte (engl. default) Werte für h und args definiert. Das bedeutet, dass die Werte h = 1e-5 und args = () verwendet werden, wenn die Funktion finite_difference ohne die Argumente h und args aufgerufen wird.

In der Funktion wird zuerst die Dimension des Punktes x0 bestimmt und in der Integer n gespeichert. Der Gradient einer reellwertigen Funktion mit $n$ Variablen ist ein Vektor der Länge $n$. Daher wird die Variable grad als ein np.ndarray der Länge n aus Nullen durch die Funktion np.zeros initialisiert.

Weil wir Gl. (1.11) auf alle $n$ Komponenten des Punktes $x^k$ anwenden müssen, ist die Nutzung einer Schleife (engl. loop) sinnvoll. Hier benutzen wir eine for-Schleife, die über alle Indizes i von 0 bis n - 1 iteriert. Beachten Sie, dass Python-Indizes bei 0 beginnen und die built-in Funktion range den Endwert n nicht einschließt.

In jeder Iteration der Schleife definieren wir zuerst den Einheitsvektor ${\hat{e}}_i$, indem wir zuerst ein Null-Array der Länge n mit np.zeros(n) erstellen und dann die $i$-te Komponente auf 1 setzen. Danach können wir den i-ten Eintrag des Gradientenarrays grad gemäß Gl. (1.11) berechnen. Der einzige Unterschied zwischen unserem Code und dieser Gleichung ist, dass wir das Zusatzargument args an die Funktion func übergeben, wobei ein Stern * vor dem Argument args steht. Die Verwendung von * wirkt an dieser Stelle als einen unären Operator, der ein Objekt in seine Bestandteile entpackt (engl. unpacking). Das bedeutet, dass die Funktion func nach dem ersten Argument die weitere Argumente im Tuple args einzelnd akzeptiert.

Nach der letzten Iteration geben wir den Gradienten in Form der Variable grad zurück.

Objektivfunktion

Als nächstes implementieren wir die Objektivfunktion, deren Wert wir minimieren wollen. Gemäß Gl. (1.8) können wir die Funktion objective_function folgendermaßen definieren:

def objective_function(
    beta: np.ndarray,
    *args: np.ndarray,
) -> float:
    concentrations: np.ndarray = args[0]
    absorbances: np.ndarray = args[1]
    return np.sum((absorbances - (beta[0] + beta[1] * concentrations))**2)

Die Signatur der Funktion objective_function folgt der Typdefinition des Arguments func in der Funktion finite_difference (Callable[[np.ndarray, Any], float]):

• Das Argument beta ist vom Typ np.ndarray,
• das Argument args ist vom Typ np.ndarray, also auch Any und
• die Funktion gibt einen Float zurück.

Hier sehen wir wieder die Verwendung von * vor dem Argument args. In diesem Fall heißt das, dass die Funktion objective_function beliebig viele weitere Argumente nach beta akzeptiert.

Die Funktion objective_function definiert zuerst die Arrays concenctrations und absorbances aus dem Argument args und gibt anschließend den Wert der Verlustfunktion der kleinsten Quadrate gemäß Gl. (1.8) zurück. Hier gibt es erneut keine große Unterschiede zwischen der mathematischen Formulierung und der programmatischen Implementierung.

Mit Hilfe der Funktionen finite_difference können wir die Funktion objective_function_gradient implementieren, die den Gradienten der Objektivfunktion berechnet:

def objective_function_gradient(
    beta: np.ndarray,
    *args: np.ndarray,
) -> np.ndarray:
    concentrations: np.ndarray = args[0]
    absorbances: np.ndarray = args[1]
    grad: np.ndarray = finite_difference(
        objective_function, 
        beta, 
        args=(concentrations, absorbances),
    )
    return grad

Die Signatur dieser Funktion ist ähnlich zu der der Funktion objective_function; lediglich der Rückgabetyp ist ein np.ndarray. Erinnern Sie sich daran, dass, da die Objektivfunktion einen $n$-dimensionalen Vektor als Argument hat, der Gradient ebenfalls ein $n$-dimensionaler Vektor ist. Nach der Definition der Arrays concenctrations und absorbances berechnen wir den Gradienten der Objektivfunktion einfach durch Aufrufen der Funktion finite_difference. Der Rückgabewert wird in der Variable grad gespeichert und zurückgegeben.

Gradientenverfahren

Anschließend implementieren wir das Gradientenverfahren. Aus Gl. (1.10) folgt, dass wir zum Starten des Verfahrens den Gradienten der Objektivfunktion func_grad, den Startpunkt x0 und die Schrittweite alpha benötigen. Als Abbruchbedingung verwenden wir eine Kombination aus der maximalen Anzahl an Iterationen und der Norm des Gradienten. Das führt zu den zusätzlichen Argumenten max_iter und max_norm. Außerdem brauchen wir noch das Argument args, das die weiteren Argumente für die Funktion func_grad enthält.

def gradient_descent(
    func_grad: Callable[[np.ndarray, Any], np.ndarray],
    x0: np.ndarray,
    alpha: float = 0.001,
    max_norm: float = 1e-6,
    max_iter: int = 10000,
    args: tuple = (),
) -> tuple[np.ndarray, int]:
    x: np.ndarray = np.copy(x0)
    for niter in range(0, max_iter):
        grad: np.ndarray = func_grad(x, *args)
        x = x - alpha * grad
        if np.linalg.norm(grad) < max_norm:
            break
    if niter == max_iter - 1:
        print('Warning: Maximum iterations reached. '
              'Result may not be reliable.')
    
    return x, niter

Zuerst kopieren wir explizit den Array x0 mit der Funktion np.copy und speichern die Kopie in der Variable x. Das ist notwendig, da Python standardmäßig Arrays nicht kopiert, sondern nur Referenzen auf sie speichert. Das bedeutet, dass, wenn wir x0 ändern, auch x geändert wird. Um das zu vermeiden, kopieren wir den Array explizit.

Danach verwenden wir eine for-Schleife, die die Variable niter von 0 bis max_iter - 1 iteriert. In jeder Iteration berechnen wir den Gradienten grad durch Aufrufen der Funktion func_grad mit den Argumenten x und args. Anschließend verwenden wir Gl. (1.10) um die Variable x zu aktualisieren. Danach berechnen wir die Norm des Gradienten mit der Funktion np.linalg.norm und prüfen, ob sie kleiner als max_norm ist. Wenn ja, brechen wir die Schleife mit dem break-Befehl ab.

Nach der letzen Iteration überprüfen wir, ob die Variable niter die maximale Anzahl an Iterationen erreicht hat. Wenn ja, bedeutet das, dass das verfahren möglicherweise nicht konvergiert ist und wir geben wir eine Warnung aus. Am Ende geben wir das Optimum x so wie die Anzahl der tatsächlich benötigten I terationen niter zurück.

Anwendung

Nun können wir den implementierten Algorithmus auf die Daten aus Kapitel 1.2 anwenden. Dazu definieren wir als erstes wieder die Arrays concentrations und absorbances:

concentrations = [
    2.125, 4.250, 6.375, 8.500, 10.63, 12.75, 14.88, 17.00, 19.13, 21.25,
    23.38, 25.50, 27.63, 29.75, 31.88, 34.00, 36.13, 38.25, 40.38, 42.50,
]
absorbances = [
    0.0572, 0.1391, 0.2049, 0.2754, 0.3420, 
    0.4139, 0.4956, 0.5815, 0.6806, 0.7481,
    0.8242, 0.9130, 1.0043, 1.0809, 1.1511,
    1.2483, 1.3373, 1.4027, 1.4927, 1.5853,
]
concentrations: np.ndarray = np.array(concentrations)
absorbances: np.ndarray = np.array(absorbances)

Dann definieren wir den Startpunkt x0, hier beta_guess, und rufen die Funktion gradient_descent auf:

beta_guess = np.array([1.0, 1.0])
beta_opt, niter = gradient_descent(
    objective_function_gradient, 
    beta_guess,
    alpha=0.00005,
    max_iter=100000,
    args=(concentrations, absorbances),
)

beta0, beta1 = beta_opt
assert np.isclose(beta0, -0.04907034)
assert np.isclose(beta1, 0.03800109)

print(beta_opt)
print(niter)

Die optimalen Parameter beta0 und beta1 sind identisch wie die der analytischen Lösung. Auf dem Computer des Autors wurden 34683 Iterationen benötigt, um die Abbruchbedingung zu erfüllen. Die genaue Anzahl der Iterationen kann je nach Hardware leicht variieren. Wählt man die Schrittweite alpha geringfügig größer, werden weniger Iterationen benötigt. Ist alpha aber zu groß, divergiert das Verfahren.

Da Optimierung ein sehr allgemeines Problem ist, existieren viele Implementierungen von verschiedensten Algorithmen in Bibliotheken wie z.B. scipy.optimize. Wir wollen die Funktion scipy.optimize.minimize verwenden, um die optimalen Parameter $\beta$ zu finden:

from scipy.optimize import minimize

res = minimize(
    objective_function,
    beta_guess,
    args=(concentrations, absorbances),
    method='CG',
    jac=objective_function_gradient,
    options={'maxiter': 10000, 'gtol': 1e-6},
)

beta0, beta1 = res.x
niter = res.nit
assert np.isclose(beta0, -0.04907034)
assert np.isclose(beta1, 0.03800109)

Beim Aufrufen der Funktion minimize müssen wir nur die Objektivfunktion objective_function und den Startpunkt beta_guess angeben. Um die Berechnung der numerischen Gradienten kümmert sich die minimize-Funktion selbst.

Mit dem Argument method='CG' haben wir die Methode des nichtlinearen Konjugierten Gradienten (engl. nonlinear Conjugate Gradient method) ausgewählt. Diese Methode ist ein Verbesseung des Gradientenverfahrens und erreicht das Minimum in nur 2 Iterationen auf dem Computer des Autors.

Bei minimize gibt es eine Reihe von weiteren Minimierungsmethoden, die verwendet werden können. Eine Übersicht finden Sie in der Dokumentation dieser Funktion. Zwei wichtige Methoden davon sind:

• method='Nelder-Mead': Das Nelder-Mead-Verfahren ist eine heutristische Methode, die ohne die Berechnung des Gradienten auskommt. Sie ist daher besonders nützlich, wenn die Berechnung des Gradienten sehr aufwendig ist oder dieser stark variiert. Sie eignet sich deshalb besonders für Regressionen mit experimentellen Daten.
• method='BFGS': Das Broyden-Fletcher-Goldfarb-Shanno-Verfahren ist eine Methode, die den Gradienten benutzt, um die Hesse-Matrix der Objektivfunktion zu approximieren. Die Hesse-Matrix enthält die zweiten Ableitungen der Funktion nach ihren Parametern. Die Methode zeigt deshaln eine sehr schnelle Konvergenz in der Nähe des Minimums. In der Praxis benögtigt sie weniger Iterationen als andere Optimierungsmethoden und wird deshalb häufig verwendet.

Wir werden in der Übung sehen, dass es für die Regression mit der Methode der kleinsten Quadrate eine geschlossene Lösung gibt, die die optimalen Parameter direkt berechnet. Warum sollten wir dann die numerische Optimierung verwenden? Im Kontext der Regression erlaubt uns die numerische Optimierung einerseits den Einsatz komplizierterer Modelle, die keine analytische Lösung haben, und andererseits die Verwendung sophistizierterer Verlustfunktionen, wie z.B. die der Methode der kleinsten absoluten Abweichungen (vgl. Gl. (1.6)). Zudem können wir damit eine zusätzliche Kontrolle über die Parameter einführen (Regularisierung), was die allgemeine Leistung des Modells verbessern kann.

Übung

Aufgabe 1.2: Polynomiale Regression

Basierend auf der vorherigen Aufgabe, in welcher Sie von der linearen Regression zur quadratischen Regression übergegangen sind, können Sie bereits vermuten, dass die Methode der kleinsten Quadrate auch mit höhergradige Polynomen simpel zu implementieren ist. In der Praxis ist es jedoch nicht sinnvoll, die Gleichungssysteme für Polynome höheren Grades manuell zu lösen. Stattdessen können Sie die Funktion np.polyfit verwenden, um die Koeffizienten ${\beta }_0,{\beta }_1,\dots ,{\beta }_n$ eines Polynoms $n$-ten Grades $$\begin{array}{c}\hat{f}(\beta ;x_i)={\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2+\dots +{\beta }_n{x}_i^n\end{array}$$ zu bestimmen, welches am besten zu den Daten passt. Diese Funktion nimmt als Argumente die Arrays der unabhängigen Variable $x$ und der abhängigen Variable $y$, sowie den Grad des Polynoms $n$ entgegen und gibt die Koeffizienten ${\beta }_j$ zurück.

(a) Polynomiale Regression mit np.polyfit

Wenden Sie die Funktion np.polyfit auf die Methylenblau-Daten an, um ein Polynom 20. Grades zu fitten und plotten Sie das Polynom zusammen mit den Datenpunkten.

(b) Vorhersage von neuen Datenpunkten

Aus dem Plot der polynomialen Regression (und ggf. den Residuen) können Sie erkennen, dass das Polynom 20. Grades die Datenpunkte sehr gut anpasst. Dies ist auch nicht weiter verwunderlich, da wir eine Funktion mit mind. 20 Parametern so anpassen können, dass sie unsere 20 Datenpunkte perfekt wiedergibt. Allerdings haben wir in unserem Code bisher lediglich die Datenpunkte, welche wir zum Fitten des Modells verwendet haben, zur Visualisierung der Ergebnisse beachtet. Die Funktionswerte zwischen den Datenpunkten wurden lediglich interpoliert. In der Regel möchten wir allerdings mit Hilfe unseres Modells auch Vorhersagen für neue Datenpunkte erhalten.

Plotten Sie das gesamte Polynom 20. Grades zusammen mit den Datenpunkten in dem Interall $0\le x\le 50.0$. Definieren Sie sich dazu ein Array mit 1000 Werten mit Hilfe der Funktion np.linspace und berechnen Sie die Funktionswerte. Beschränken Sie die Darstellung des Plots auf den Bereich $0\le y\le 2.0$. Was beobachten Sie?

Nichtlineare Regression

Wir haben in den vorherigen Abschnitten die lineare Regression kennengelernt, die uns erlaubt, lineare Zusammenhänge zwischen Variablen zu modellieren. Zwar sind viele physikalische Zusammenhänge linear, oder können als solche formuliert werden, aber es gibt auch viele nichtlineare Zusammenhänge, die wir modellieren wollen. Im Rahmen der Methode der kleinsten Quadrate ersetzen wir dazu einfach das Modell $\widehat{f(\beta ;x)}$ in Gl. (1.4) durch eine nichtlineare Funktion. Allerdings ist in diesem Fall ist eine analytische Lösung wie Gl. (1.9) nicht immer möglich, weshalb numerische Optimierungsverfahren verwendet werden müssen.

Anwendung

Reaktionskinetik

Sie haben im Physikalisch-Chemischen Praktikum sicherlich den Versuch “Bestimmung der Geschwindigkeitskonstante und der Aktivierungsenergie der Mangan(III)-Trioxalat-Zersetzungsreaktion”, auch “Mn-Zerfall” genannt, durchgeführt. Dort haben Sie die Absorbanz $A$ in Abhängigkeit der Zeit $t$ gemessen und durch fitten der Messdaten die Geschwindigkeitskonstante $k$ bestimmt. Die zugrungeliegende Beziehung ist exponentiell: $$A(t)=A_0{\mathrm{e}}^{-kt}$$ mit dem Parametern $A_0$ und $k$, d.h. $\beta =(A_0,k{)}^{⊺}$.

Hier gilt also $A(t)=\hat{f}(\beta ;t)$ und wir können mit Hilfe der Verlustfunktion der kleinsten Quadrate das Regressionsproblem als das folgende Optimierungsproblem formulieren: $${\beta }^{\ast }=\underset{\beta \in {\mathbb{R}}^2}{\mathrm{argmin}}\sum _{i=1}^N(A_i-A_0{\mathrm{e}}^{-k{t}_i}{)}^2\text{\hspace{1em}(1.12)}$$

Wir importieren als erstes wieder die benötigten Module und Bibliotheken:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

Anschließend müssen wir die Daten in Form von Arrays bereitstellen. Weil wir hier doch relativ viele Datenpunkte haben, wird das manuelle Eintippen ziemlich mühsam. Deshalb verwenden wir die Funktion np.loadtxt, um die Daten aus einer Textdatei zu lesen. Die Textdatei mn_decay.txt (hier herunterladen) enthält zwei Spalten, welche die Werte für Zeit $t$ und Absorbanz $A$ enthalten. Die ersten Zeilen der Datei sehen wie folgt aus:

# time / s  absorbance
         0    1.625020
        30    1.471490
        60    1.295880
        90    1.169340
       120    1.042040
       150    0.924193
       180    0.824484
       210    0.731743
       240    0.650593

Zum Einlesen der Daten benötigt die Funktion np.loadtxt den Dateinamen:

time, absorbance = np.loadtxt('mn_decay.txt', unpack=True)

In diesem Fall liegt die Textdatei im gleichen Verzeichnis wie das ausführende Skript. Wenn Sie die Textdatei in einem anderen Verzeichnis haben, müssen Sie den Pfad entsprechend anpassen. Das optionale Argument unpack=True sorgt dafür, dass die Spalten der Datei einzeln ausgegeben und als Arrays time und absorbance gespeichert werden. Würden wir unpack=False setzen, welches auch der Defaultwert ist, würde die Funktion ein 2D-Array zurückgeben, in dem die Spalten zusammengefasst sind. Die erste Zeile dieser Datei beginnt mit einem Kommentarzeichen #, was np.loadtxt dazu veranlasst, diese Zeile zu ignorieren.

Nun können wir das Modell definieren:

def exp_decay(t: np.ndarray, a0: float, k: float) -> np.ndarray:
    return a0 * np.exp(-k * t)

Obwohl t ein Array und k ein Skalar ist, funktioniert die Multiplikation -k * t elementweise. In diesem Fall wird auch k als ein Array interpretiert, was als Broadcasting bezeichnet wird. Die Funktion np.exp berechnet den elementweisen Exponential des Arrays. Das Ergebnis ist demnach wieder ein Array, welches wir mit dem Skalar a0 multiplizieren.

Anschließend definieren wir die Objektivfunktion in Gl. (1.12):

def objective_function(beta, *args):
    a0, k = beta
    time, absorbance = args
    return np.sum((absorbance - exp_decay(time, a0, k))**2)

Wir verwenden nun die minimize-Funktion, um dieses Optimierungsproblem zu lösen:

beta_guess = (1.0, 0.01)
res = minimize(
    objective_function, beta_guess, method='Nelder-Mead', 
    args=(time, absorbance),
)
a0, k = res.x

print(f'k = {k} s^-1')
print(f'A0 = {a0}')

Hier haben wir das Nelder-Mead-Verfahren mit den Startparametern $A_0^0=1$ und $k^0=0.01$ angewandt. Zum Ausgeben des Ergebnisses mit dem print-Befehl haben wir vor der Zeichenkette jeweils ein f gesetzt. Das signalisiert, dass die Zeichenkette ein sog. f-string ist, in welchen wir Variablen mit geschweiften Klammern {} einbetten können. Tatsächlich können f-Strings noch einiges mehr, was wir in Zukunft noch sehen werden. Die optimierten Parameter sollten die folgenden Werte haben:

assert np.isclose(k, 0.0039022970)
assert np.isclose(a0, 1.6475263)

Zum Schluss können wir die Ergebnisse plotten:

fig, ax = plt.subplots(figsize=(8, 6))
time_interp = np.linspace(time.min(), time.max(), 1000)

ax.plot(time, absorbance, 'o', label='data')
ax.plot(time_interp, exp_decay(time_interp, a0, k), label='exponential fit')

ax.set_xlabel('time / s')
ax.set_ylabel('absorbance')
ax.legend()

fig.tight_layout()
plt.show()

Sie sollten die meisten Funktionen im obigen Codeblock aus Kap. 1.2 kennen. Ein Unterschied ist die Verwendung der Funktion np.linspace, welche Zeiten zwischen den Messpunkten generiert, sodass wir das Regressionsmodell für eine Interpolation verwenden können. Diese Funktion akzeptiert drei Argumente: den Startwert, den Endwert und die Anzahl der zu generierenden Punkte. Dann produziert sie ein Array mit gleichmäßig verteilten Werten zwischen dem Start- und Endwert.

Ein weitere neue Funktion ist fig.tight_layout(), die eine automatische Anpassung des Layouts des Plots vornimmt. Wir erkennen aus dem folgendem Diagramm, dass die exponentielle Funktion die Daten sehr gut beschreibt.

Einige von Ihnen würden vielleicht fragen, warum wir nicht die Funktion linearisiert haben, um die lineare Regression zu verwenden, was eine berechtigte Frage ist. Tatsächlich ist es möglich, die Funktion zu linearisieren, indem wir beide Seiten der Gleichung logarithmieren: $$\ln (A(t))=\ln (A_0)-kt$$ Eine lineare Regression mittels dieser Gleichung liefert allerdings nicht die gleichen Ergebnisse, wie Sie im folgenden Diagramm sehen können:

Als freiwillige Übung können Sie versuchen, das obige Diagramm reproduzieren. Es ist unschwer zu erkennen, dass der linearisierte Fit schlechter zu den Daten passt. Das liegt daran, dass die lineare Regression die Fehler in der Absorbanz durch das Logarithmieren nicht gleichmäßig behandelt. Die dadurch erhaltenen Parameter

assert np.isclose(k_lin, 0.0041675912)
assert np.isclose(a0_lin, 1.7633374)

sind durchaus unterschiedlich zu den vorherigen. Deshalb ist es oft notwendig, nichtlineare Regressionen an den urprünglichen Daten durchzuführen, anstatt lineare Modelle mit linearisierten Daten zu verwenden.

Titrationskurve

Im Analytikpraktikum haben Sie sicherlich ebenfalls eine Titration einer starken Base gegen eine starke Säure mit einem pH-Meter durchgeführt. Damals mussten Sie die Werte wahrscheinlich auf einem Millimeterpapier auftragen und anhand der Position des pH-Sprungs den Äquivalenzpunkt bestimmen. Das ist einerseits mühsam und andererseits ungenau, da nur die wenigen Messdaten in der Nähe des steilen Anstiegs berücksichtigt werden.

Da die pH-Kurve eine Funktion in Abhängigkeit der zugegebenen Menge an Base ist, können wir sie mithilfe der nichtlinearen Regression modellieren und den Äquivalenzpunkt mit deutlich höherer Genauigkeit bestimmen.

Die ${\mathrm{H}}^{+}$-Konzentration während der Titration einer starken Base gegen eine starke Säure ist (unter gewissen Näherungen) gegeben durch: $$[{\mathrm{H}}^{+}]=\frac{\Delta +\sqrt{{\Delta }^2+4K_w}}{2},\text{\hspace{1em}(1.13)}$$ wobei $\Delta =[{\mathrm{A}}^{-}]-[{\mathrm{B}}^{+}]$ die Konzentrationsdifferenz zwischen den Gegenionen der Säure ${\mathrm{A}}^{-}$ und der Base ${\mathrm{B}}^{+}$ ist. $K_w$ ist das Ionenprodukt des Wassers.

Da starke Säuren und Basen vollständig dissoziieren, lassen sich die Konzentrationen ihrer Gegenionen wie folgt ausdrücken: $$\begin{array}{rl}[{\mathrm{A}}^{-}]& =\frac{n_{\mathrm{A}}}{V}=\frac{c_{\mathrm{A}}^0{V}^0}{V^0+V_{\mathrm{B}}}\\ [{\mathrm{B}}^{+}]& =\frac{n_{\mathrm{B}}}{V}=\frac{c_{\mathrm{B}}^0{V}_{\mathrm{B}}}{V^0+V_{\mathrm{B}}},\end{array}$$ wobei $c_{\mathrm{A}}^0$ und $c_{\mathrm{B}}^0$ die Konzentrationen der zu analysierenden Säure und der zugegebenen Base sind, $V^0$ das Anfangsvolumen der Probelösung und $V_{\mathrm{B}}$ das Volumen der zugegebenen Base ist.

Der pH-Wert lässt sich aus der ${\mathrm{H}}^{+}$-Konzentration berechnen: $$\mathrm{pH}=-\lg \left(\frac{[{\mathrm{H}}^{+}]}{1\mathrm{M}}\right).\text{\hspace{1em}(1.14)}$$

Fasst man die Gleichungen (1.13) und (1.14) in einer Funktion zusammenfassen, so erhalten wir das Modell $f(\beta ;V_{\mathrm{B}})$, wobei $\beta =(c_{\mathrm{A}}^0,V^0{)}^{⊺}$.

Wir implmentieren zunächst das Modell und die Objektivfunktion:

def titration_sasb_model(
    v_b: np.ndarray,
    c0_b: float,
    c0_a: float,
    v0: float,
) -> np.ndarray:
    k_w = 1e-14
    c_a = c0_a * v0 / (v0 + v_b)
    c_b = c0_b * v_b / (v0 + v_b)
    delta = c_a - c_b
    c_h = 0.5 * (delta + np.sqrt(delta**2 + 4 * k_w))
    ph = -np.log10(c_h)
    return ph

def objective_function(beta, *args):
    c0_a, v0 = beta
    c0_b, v_b, ph = args
    ph_fit = titration_sasb_model(v_b, c0_b, c0_a, v0)
    return np.sum((ph - ph_fit)**2)

Obwohl die Funktion des pH-Werts relativ kompliziert ist, können wir durch die Definition von Zwischenvariablen, wie in Gl. (1.13), die Implementierung der Funktion in Python deutlich vereinfachen. Die Objektivfunktion ist fast identisch zu der des Mn-Zerfalls, wobei der wesentliche Unterschied die Ersetzung unseres Modells darstellt.

Genau so wie im vorherigen Beispiel lesen wir die Daten aus einer Textdatei (hier herunterladen) ein:

C0_B = 0.1  # mol/l
v_b, ph = np.loadtxt('titration_sasb.txt', unpack=True)

Zusätzlich haben wir hier die Konzentration der Maßlösung C0_B definiert. Gemäß der allgemeinen Konvention sollen alle Konstanten in Python in Großbuchstaben geschrieben werden. Anschließend können wir die nichtlineare Regression durchführen und die Ergebnisse plotten:

beta_guess = (0.01, 100)
res = minimize(
    objective_function, 
    beta_guess, 
    args=(C0_B, v_b, ph),
    method='Nelder-Mead',
)

c0_a, v0 = res.x
print(f'c0_a = {c0_a} mol/l')
print(f'v0 = {v0} ml')
print(f'n0_a = {c0_a * v0} mmol')

fig, ax = plt.subplots(figsize=(8, 6))
v_b_interp = np.linspace(v_b.min(), v_b.max(), 1000)

ax.plot(v_b, ph, 'o', label='data')
ax.plot(
    v_b_interp, 
    titration_sasb_model(v_b_interp, C0_B, c0_a, v0), 
    label='pH fit',
)

ax.set_xlabel('volume of base / ml')
ax.set_ylabel('pH')
ax.legend()

fig.tight_layout()
plt.show()

Aus den gefitteten Parametern wurde eine Stoffmenge des Analyts von $n_{\mathrm{A}}^0=1.7698\mathrm{mmol}$ bestimmt. Das entsprechende Diagramm sollte wie folgt aussehen: Obwohl der Fit am Anfang und am Ende der Kurve nicht perfekt ist, ist die Übereinstimmung in der Nähe des Äquivalenzpunkts sehr gut.

Vielleicht haben Sie sich gefragt, warum die optimierten Parameter $c_{\mathrm{A}}^0$ und $V^0$ hier nicht explizit aufgeführt sind, sondern lediglich ihr Produkt.

Als letztes betonen wir nochmal, dass es oft wichtig ist, die Regression an den Originaldaten durchzuführen und keine transformierten Daten zu verwenden. Im Fall der Säure-Base-Titration könnte man z.B. auf die Idee kommen, $[{\mathrm{H}}^{+}]$ anstatt des pH-Werts zu fitten. Aus den gleichen Gründen wie zuvor ist das allerdings nicht sinnvoll: Da der gleiche Fehler auf der pH-Skala zu größeren Fehlern bei höheren ${\mathrm{H}}^{+}$-Konzentrationen und kleineren Fehlern bei niedrigeren ${\mathrm{H}}^{+}$-Konzentrationen führt, werden die Datenpunkte bei höheren ${\mathrm{H}}^{+}$-Konzentrationen besser angepasst. Das führt zu einer Verzerrung der Ergebnisse, wie im folgenden Diagramm zu sehen ist: Als freiwillige Übung können Sie versuchen, das obige Diagramm reproduzieren. Man erkennt, dass die Regression an $[{\mathrm{H}}^{+}]$ die früheren Datenpunkte bevorzugt und dadurch den Äquivalenzpunkt völlig falsch bestimmt wird.

Übung

Aufgabe 1.3: Regularisierung

Das Phänomen, welches Sie bei der polynomialen Regression 20. Ordnung beobachten können, wird als Überanpassung (engl. overfitting) bezeichnet. Es tritt auf, wenn das Modell zu komplex ist und nicht nur der zugrunde liegende Trend, sondern auch das Rauschen in den Daten angepasst wird. In solchen Fällen kann das Modell die Datenpunkte zwar perfekt reproduzieren, aber es wird nicht in der Lage sein, neue Datenpunkte vorherzusagen.

Um Überanpassung zu vermeiden gibt es, neben der Reduzierung der Parameter, die Möglichkeit der Regularisierung. Darunter versteht man die Einführung von zusätzlichen Bedingungen, welche die Komplexität des Modells einschränken. Eine solche Bedingung kann beispielsweise sein, dass die Koeffizienten ${\beta }_j$ möglichst klein gehalten werden, was durch die Einführung eines zusätzlichen Terms in die Verlustfunktion erreicht werden kann.

Verwendet man das Quadrat der ${\ell }_2$-Norm der Koeffizienten $\Vert \beta {\Vert }^2={\sum }_i{\beta }_i^2$ als Regularisierung und fügt sie der Verlustfunktion hinzu, so spricht man von Ridge-Regression. Die Verlustfunktion ist dann gegeben durch $$\begin{array}{c}L(\beta ;x,y)=\sum _{i=1}^N(y_i-\hat{f}(\beta ;x_i){)}^2+\lambda \Vert \beta {\Vert }^2,\end{array}$$ wobei der Parameter $\lambda$ die relative Stärke der Regularisierung bestimmt.

(a) Ridge-Regression der Methylenblau-Daten mit Polynom 20. Ordnung

Implementieren Sie die Ridge-Regression für die Methylenblau-Daten mit $\lambda =0.001$ und fitten Sie ein Polynom 20. Ordnung. Nutzen Sie dazu die numerische Optimierungsmethode mit der Funktion minimize und ändern Sie Ihre Objektivfunktion entsprechend. Verwenden Sie als Startwerte ein Array mit Nullen. Normalisieren Sie außerdem vor der Regression die Konzentrationen und die Absorptionswerte auf den Bereich $[0,1]$, indem Sie jeweils durch den Maximalwert teilen. Plotten Sie das Ergebnis zusammen mit den Datenpunkten.

(b) Einfluss des Regularisierungsparameters $\lambda$

Variieren Sie den Regularisierungsparameter $\lambda$ und beobachten Sie, wie sich die Stärke der Regularisierung auf die Anpassung des Modells an die Datenpunkte auswirkt. Was passiert, wenn Sie $\lambda =0$ oder $\lambda =1$ wählen?

Differentialgleichungen

Im letzten Kapitel haben wir mit Hilfe der nichtlinearen Regression den zeitlichen Verlauf der Konzentration eines Mn-Komplexes untersucht. Dies ist ein Musterbeispiel für eine Reaktion erster Ordnung, welche durch eine vergleichsweise einfache zeitliche Entwicklung charakterisiert ist. Allerdings gibt es viele weitere chemische Reaktionen, bei denen die Dynamik deutlich komplexer ist. Zu diesem Zweck benötigen wir mathematische Modelle, die die zeitliche Entwicklung von chemischen Systemen beschreiben können, und welche auf Differentialgleichungen (DGLs) basieren.

Ein weiteres wichtiges Anwendungsgebiet von DGLs in der Chemie ist die Beschreibung quantenchemischer Systeme durch die Schrödingergleichung. Sie wissen sicherlich, dass die Schrödingergleichung für Modellsysteme wie das Teilchen im Kasten oder der harmonischen Oszillator analytisch lösbar ist. Für molekulare Systeme ist dies jedoch nicht der Fall. Die Quantenchemie basiert daher auf der (approximativen) numerischen Lösung dieser Gleichung, um die elektronischen Struktur und Eigenschaften von Molekülen zu berechnen.

In diesem Kapitel werden wir uns auf numerische Methoden zur Lösung von Differentialgleichungen konzentrieren, da diese in der Chemie eine wichtige Rolle spielen. Wir werden verschiedene numerische Verfahren kennenlernen und diese auf praktische Beispiele aus der chemischen Kinetik und Quantenchemie anwenden. Durch die Anwendung dieser Methoden werden Sie in die Lage versetzt, komplexe chemische Systeme zu modellieren und die Dynamik dieser Systeme zu verstehen.

Anfangswertproblem

Eine Differentialgleichung ist eine Gleichung, die eine unbekannte Funktion und ihre Ableitungen enthält. Eine allgemeine Form solcher Gleichungen ist kompliziert, da die Funktion und ihre (partiellen) Ableitungen von verschiedester Ordnung und zu unterschiedlichen Potenzen auftreten können. In diesem Abschnitt werden wir uns auf die gewöhnlichen Differentialgleichungen (GDGLs, engl. ODEs) beschränken, bei denen die Funktion nur von einer Variablen abhängt. Außerdem werden wir nur lineare Differentialgleichungen betrachten, bei denen die Funktion und ihre Ableitungen nur in der ersten Potenz auftreten (d.h. ersten Grades sind). Die allgemeine Form einer linearen gewöhnlichen Differentialgleichung lautet: $$a_n(x)y^{(n)}(x)+a_{n-1}(x)y^{(n-1)}(x)+\dots +a_1(x)y^{\prime }(x)+a_0(x)y(x)=b_0(x)u(x),\text{\hspace{1em}(2.1)}$$ wobei $x$ die unabhängige Variable, $y(x)$ die gesuchte Funktion, $y^{(i)}(x)$ die $i$-te Ableitung dieser Funktion ist und $a_i(x)$, $b_0(x)$ und $u(x)$ gegebene Funktionen von $x$ sind. Das größte $i$, für das $a_i(x)\ne 0$, ist die Ordnung der Differentialgleichung. Eine lineare Differentialgleichung erster Ordnung hat demnach die Form $$a_1(x)y^{\prime }(x)+a_0(x)y(x)=b_0(x)u(x).\text{\hspace{1em}(2.2)}$$

Eine allgemeine DGL erster Ordnung kann in der expliziten Form $$y^{\prime }(x)=f(x,y(x))\text{\hspace{1em}(2.3)}$$ geschrieben werden, wobei $f(x,y(x))$ eine gegebene Funktion von der unabhängigen Variablen $x$ und der gesuchten Funktion $y(x)$ ist. Es ist offensichtlich, dass Gl. (2.3) ein spezieller Fall von Gl. (2.2) ist, mit $f(x,y(x))=-\frac{a_0(x)}{a_1(x)}y(x)+\frac{b_0(x)}{a_1(x)}u(x)$.

Da bei der Bildung der Ableitung der konstante Term verschwindet, ist die (allgemeine) Lösung $y(x)$ einer DGL nicht eine eindeutige Funktion, sondern eine Funktionsschar, die durch Integrationskonstanten parametrisiert ist. Um eine eindeutige Lösung zu erhalten, müssen diese Konstanten bestimmt werden, wie z.B. durch die Wahl von Anfangsbedingungen.

Ein Anfangswertproblem (AWP) ist dann ein Problem, bei dem zusätzlich zu der Differentialgleichung eine oder mehrere Anfangsbedingungen gegeben sind. Für eine lineare GDL erster Ordnung ist die Anfangsbedingung durch den Funktionswert $y(x_0)=y_0$ gegeben, wobei $x_0$ ein gegebener Punkt und $y_0$ ein gegebener Funktionswert sind. Obwohl wie dies als Anfangsbedingung bezeichnen, muss der gegebene Punkt nicht notwendigerweise der Anfangspunkt, z.B. Zeitpunkt, sein. Die Lösung eines AWP ist dann eine eindeutige Funktion, die auch als spezielle Lösung der DGL bezeichnet wird.

Für eine DGL $n$-ter Ordnung benötigen wir $n$ Anfangsbedingungen, um eine spezielle Lösung zu erhalten. Diese sind gegeben durch $$y(x_0)=y_0,y^{\prime }(x_0)=y_0^{\prime },\dots ,y^{(n-1)}(x_0)=y_0^{(n-1)}.$$

Es gibt noch weitere Möglichkeiten, eine spezielle Lösung zu erhalten, z.B. durch die Angabe von Randbedingungen, was Sie in der Übung kennenlernen werden.

Euler-Verfahren

Wir betrachten das Anfangswertproblem $$\{\begin{array}{l}{y}^{\prime }(x)=f(x,y(x))\\ y(x_0)=y_0\end{array},\text{\hspace{1em}(2.4)}$$ wobei $f(x,y(x))$ eine gegebene Funktion von der unabhängigen Variablen $x$ und der gesuchten Funktion $y(x)$ ist. Gl. (2.4) ist ein AWP erster Ordnung, mit einer DGL erster Ordnung wie in Gl. (2.3) und einer Anfangsbedingung.

Unser Ziel ist es, die spezielle Lösung $y(x)$ des AWP numerisch zu finden. Wie bei vielen numerischen Verfahren, beginnen wir mit einer Diskretisierung der Funktion $y(x)$, d.h. wir wählen eine Menge von Punkten $x_i$ und betrachten die Funktion nur an diesen Punkten anstatt auf dem gesamten Definitionsbereich. Die konzeptionell wohl einfachste Wahl der Punkte ist ein gleichmäßiges Gitter (Grid), was bedeutet, dass man einen Anfangspunkt $x_0$ und eine Schrittweite $h$ wählt, so dass die weiteren Punkte durch $x_n=x_0+nh$ für $n=0,1,2,\dots$ festgelegt werden. Unser Ziel in den folgenden Abschnitten wird es sein, die Funktionswerte von $y(x)$ an diesen Punkten $x_1,x_2,\dots$ zu berechnen, bzw. zu approximieren.

Der Funktionswert von $y$ an dem Punkt $x_{n+1}$ kann (unter bestimmten Voraussetzung an $y(x)$) durch eine Taylor-Entwicklung um den Punkt $x_n$ ersetzt werden, also $$y(x_{n+1})=y(x_n)+h{y}^{\prime }(x_n)+\frac{h^2}{2}{y}^{\prime \prime }(x_n)+\frac{h^3}{6}{y}^{\prime \prime \prime }(x_n)+\cdots .$$

Theoretische Grundlagen

Gehen wir nun davon aus, dass die Funktion $y(x)$ gut durch ihr Taylor-Polynom 1. Ordnung $$y(x_{n+1})=y(x_n)+h{y}^{\prime }(x_n)+\mathcal{O}(h^2)$$ approximiert werden kann, also dass der Fehler $\mathcal{O}(h^2)$, der proportional zu $h^2$ ist, klein ist. Mathematisch unsauber, aber praktisch für die Implementierung, schreiben wir im Folgenden $$y(x_{n+1})=y(x_n)+h{y}^{\prime }(x_n),\text{\hspace{1em}(2.5)}$$ wobei wir den Fehler vernachlässigen. Gl. (2.5) zeigt, dass wir, wenn wir den Funktionswert $y$ und die Ableitung $y^{\prime }$ an dem Punkt $x_n$ kennen, den Funktionswert an dem nächsten Punkt $x_{n+1}$ berechnen können. Mit Hilfe der DGL in Gl. (2.4) können wir dabei die Ableitung $y^{\prime }(x_n)$ durch $f(x_n,y(x_n))$ ersetzen, was zu $$y(x_{n+1})=y(x_n)+hf(x_n,y(x_n))$$ führt. Um zu betonen, dass die Euler-Methode nur diskrete Werte von $y(x)$ liefert, schreiben wir $y_n=y(x_n)$ und $y_{n+1}=y(x_{n+1})$, was zu $$y_{n+1}=y_n+hf(x_n,y_n)\text{\hspace{1em}(2.6)}$$ führt. Gl. (2.6) beschreibt das explizite Euler-Verfahren für die numerische Lösung eines AWP erster Ordnung. Durch die Anfangsbedingung kennen wir $y(x_0)$ und können dann, Schritt für Schritt, alle weiteren $y(x_n)$ mit Hilfe von Gl. (2.6) berechnen.

Implementierung

Mn-Zerfall

Wir nehmen wieder den Mn-Zerfall als Beispiel. Dort haben wir eine Reaktion 1. Ordnung, welche durch die DGL $$c^{\prime }(t)=-kc(t)\text{\hspace{1em}(2.7)}$$ beschrieben wird, wobei $c(t)$ die Konzentration des Mn-Komplexes zum Zeitpunk $t$ beschreibt. Wir können nach dem Importieren der benötigten Libraries

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable

die Funktion dydx implementieren, die $f(x_n,y(x_n))$ berechnet:

def dydx(x: float, y: float) -> float:
    k = 0.0039022970
    return -k * y

Diese Funktion akzeptiert die Argumente x ($t_n$) und y ($c(t_n)$) und gibt die Ableitung $c^{\prime }(t_n)=-kc(t_n)$ zurück. Hier haben wir die gefittete Geschwindigkeitskonstante $k$ aus Abschnitt 1.4 verwendet. Anschließend können wir gemäß Gl. (2.6) die Funktion euler_step implementieren, die den Funktionswert $y(x_{n+1})$ berechnet:

def euler_step(
    x_n: float,
    y_n: float,
    h: float,
    dydx: Callable[[float, float], float],
) -> float:
    return y_n + h * dydx(x_n, y_n)

Wir können jetzt das Euler-Verfahren implementieren:

def euler_method(
    x0: float, 
    y0: float, 
    h: float, 
    dydx: Callable[[float, float], float], 
    nsteps: int,
) -> np.ndarray:
    x = x0 + np.arange(0, nsteps + 1) * h
    y = np.zeros(nsteps + 1)
    y[0] = y0

    for i in range(0, nsteps):
        y[i + 1] = euler_step(x[i], y[i], h, dydx)

    return x, y

Diese Funktion akzeptiert neben der Anfangsbedingungen x0 und y0 die Schrittweite h, die Ableitungsfunktion dydx und die Anzahl der Schritte n. Wir erstellen zunächst das Grid x mit der Funktion np.arange und initialisieren das Nullarray y, um später die Lösung zu speichern. Dann wird der erste Eintrag y[0] dieses Arrays mit dem Anfangswert y0 überschrieben. Anschließend verwenden wir eine for-Schleife über die Anzahl der Schritte und rufen in jedem Schritt die Funktion euler_step auf, die den Funktionswert an dem jeweils nächsten Punkt berechnet, und speichern diesen in y[i + 1]. Am Ende wird das Grid x und die Lösung y zurückgegeben.

Nun wenden wir das Euler-Verfahren auf Gl. (2.7) an:

C0 = 1.0  # M
T0 = 0.0
STEP = 1.0  # s
MAXTIME = 900.0  # s

nsteps = int(MAXTIME / STEP)
x, y = euler_method(T0, C0, STEP, dydx, nsteps)

Wir setzen dazu zunächst die Anfangsbedingungen C0 = 1.0 und T0 = 0.0, die Schrittweite h = 1.0 sowie die maximale Zeit MAXTIME = 600.0. Die Anzahl der Schritte nsteps wird durch int(MAXTIME / h) berechnet. Die int-Funktion rundet dabei das Ergebnis der Division ab und konvertiert es in eine Ganzzahl. Anschließend rufen wir die Funktion euler_method auf und speichern das Ergebnis in x und y.

Zum Schluss können wir das numerische Ergebnis mit der analytischen Lösung vergleichen. Dabei sei angemerkt, dass die analytische Lösung von Gl. (2.7) leicht zu berechnen ist und durch $c(t)=c_0\exp (-kt)$ gegeben ist.

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(x, C0 * np.exp(-0.0039022970 * x), label='analytical solution')
ax.plot(x, y, label='numerical solution')

axins = ax.inset_axes(
    [0.4, 0.5, 0.27, 0.47],                 # position and size of the inset
    xlim=(300, 305), ylim=(0.3, 0.315),     # limits of the inset
    xticks=[], yticks=[]                    # remove ticks
    )
axins.plot(x, C0 * np.exp(-0.0039022970 * x))
axins.plot(x, y)
ax.indicate_inset_zoom(axins, edgecolor="black")

ax.set_xlabel('time / s')
ax.set_ylabel('concentration / M')
ax.legend(loc='upper right')

fig.tight_layout()
plt.show()

Dabei sollte dieses Diagramm erscheinen:

Erst durch Vergrößerung des Konzentrationsverlaufes können wir den Unterschied zwischen der analytischen und der numerischen Lösung erkennen, was hier mit Hilfe des Befehls ax.inset_axes erreicht wurde (die genaue Funktionsweise dieses Befehls ist an dieser Stelle unwichtig). Das Euler-Verfahren mit $h=1$ liefert demnach eine sehr gute Approximation.

Belousov-Zhabotinsky-Reaktion

Nun wollen wir eine Reaktion mit deutlich komplizierterer Kinetik betrachten. Die Belousov-Zhabotinsky-Reaktion ist ein klassisches Beispiel einer oszillierenden Reaktion, wobei ein Redox-System ($\mathrm{C}{\mathrm{e}}^{3+}/C{\mathrm{e}}^{4+}$ in der Originalreaktion) abwechselnd in der oxidierten und der reduzierten Form vorliegt. In dem folgenden Video können Sie die Reaktion beobachten:

Hier wird allerdings als Redox-System Ferroin verwendet, welches eine stärkere Farbänderung zeigt.

Der detaillierte Mechanismus dieser Reaktion ist sehr kompliziert, weshalb wir hier nur eine vereinfachte Version, das sog. Oregonator-Modell, betrachten. Ein häufig verwendetes Oregonator-Modell besteht aus fünf gekoppelten Reaktionen mit sechs Spezies:^{[1, 2]}

$$\begin{array}{rl}\mathrm{A}+\mathrm{Y}& \stackrel{k_1}{\to }\mathrm{X}+\mathrm{P}\\ \mathrm{X}+\mathrm{Y}& \stackrel{k_2}{\to }2\mathrm{P}\\ \mathrm{A}+\mathrm{X}& \stackrel{k_3}{\to }2\mathrm{X}+\mathrm{Z}\\ 2\mathrm{X}& \stackrel{k_4}{\to }\mathrm{A}+\mathrm{P}\\ \mathrm{B}+\mathrm{Z}& \stackrel{k_5}{\to }\mathrm{Y},\end{array}$$ wobei $\mathrm{A}=\mathrm{Br}{\mathrm{O}}_3^{-}$, $\mathrm{B}=\mathrm{Brommalons}\ddot{\mathrm{a}}\mathrm{ure}$, $\mathrm{X}=\mathrm{HBr}{\mathrm{O}}_2$, $\mathrm{Y}=\mathrm{B}{\mathrm{r}}^{-}$, $\mathrm{Z}=2\mathrm{C}{\mathrm{e}}^{4+}$ und $\mathrm{P}=\mathrm{HOBr}$. Möglicherweise fällt Ihnen auf, dass die obigen Reaktionsgleichungen nicht ausbalanciert sind. Der Grund dafür ist, dass die Konzentrationen einiger Spezies entweder aufgrund ihrer hohen Konzentration (${\mathrm{H}}^{+}$, Malonsäure, etc.) oder ihrer schnellen Reaktionen ($\mathrm{Br}{\mathrm{O}}_2^{\cdot }$) als konstant angenommen werden können. Insbesondere wird angenommen, dass die Konzentrationen $[\mathrm{A}]$ und $[\mathrm{B}]$ konstant sind.

Das führt zu einem System von drei gekoppelten DGLs: $$\begin{array}{rl}[X{]}^{\prime }(t)& =k_1{c}_A^0[Y]-k_2[X][Y]+k_3{c}_A^0[X]-2k_4[X{]}^2\\ [Y{]}^{\prime }(t)& =-k_1{c}_A^0[Y]-k_2[X][Y]+k_5{c}_B^0[Z]\\ [Z{]}^{\prime }(t)& =k_3{c}_A^0[X]-k_5{c}_B^0[Z]\end{array}\text{\hspace{1em}(2.8)}$$

Ein wichtiger Unterschied zum Mn-Zerfall ist, dass wir hier ein System von DGLs vorliegen haben. Glücklicherweise ist Gl. (2.6) für DGL-Systeme genauso gültig wie für einzelne DGLs, sofern wir $y(x)$ als eine vektorwertige Funktion ansetzten, d.h. $y(x)\widehat{=}([X](t),[Y](t),[Z](t){)}^{⊺}$.

Nachdem wir die benötigten Libraries importiert haben

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable

können wir die Funktion dydx für das Oregonator-Modell implementieren:

def dydx(x: float, y: np.ndarray) -> np.ndarray:
    # concentrations adapted from 
    # R. J. Field, H.-D. Försterling, J. Phys. Chem. 1986, 90, 5400–5407.
    k1 = 1.3    # M^-1 s^-1
    k2 = 2.4e6  # M^-1 s^-1
    k3 = 34.0   # M^-1 s^-1
    k4 = 3.0e3  # M^-1 s^-1
    k5 = 1.0    # M^-1 s^-1
    c_a = 0.1   # M
    c_b = 0.4   # M
    
    c_x, c_y, c_z = y

    dcxdt = k1 * c_a * c_y - k2 * c_x * c_y + k3 * c_a * c_x - 2.0 * k4 * c_x**2
    dcydt = -k1 * c_a * c_y - k2 * c_x * c_y + k5 * c_b * c_z
    dczdt = k3 * c_a * c_x - k5 * c_b * c_z

    return np.array([dcxdt, dcydt, dczdt])

Beachten Sie, dass der Datentyp des Arguments y ein np.ndarray ist. Die Funktionen euler_step und euler_method sind analog zu den Funktionen für den Mn-Zerfall implementiert, nur dass wir hier die Typ-Deklarationen in den Signaturen anpassen müssen.

def euler_step(
    x_n: float,
    y_n: np.ndarray,
    h: float,
    dydx: Callable[[float, np.ndarray], np.ndarray],
) -> np.ndarray:
    return y_n + h * dydx(x_n, y_n)

def euler_method(
    x0: float, 
    y0: np.ndarray,
    h: float, 
    dydx: Callable[[float, np.ndarray], np.ndarray],
    nsteps: int,
) -> np.ndarray:
    ndim = len(y0)

    x = x0 + np.arange(0, nsteps + 1) * h
    y = np.zeros((ndim, nsteps + 1))
    y[:, 0] = y0

    for i in range(0, nsteps):
        y[:, i + 1] = euler_step(x[i], y[:, i], h, dydx)

    return x, y

Ein weiterer Unterschied ist, dass die Variable y in der Funktion euler_method als ein 2D-Array initialisiert wird, wobei die erste Dimension die Anzahl der Komponenten angibt (hier: ndim = 3) und die zweite Dimension die Anzahl der Schritte enthält. Das Grid x bleibt ein 1D-Array, da die Zeit für alle Komponenten gleichermaßen gilt.

Nun lösen wir das DGL-System in Gl. (2.8) mit dem Euler-Verfahren:

CX_0 = 0.0    # M
CY_0 = 0.001  # M
CZ_0 = 0.0    # M
C0 = np.array([CX_0, CY_0, CZ_0])

T0 = 0.0
STEP = 1.0
TMAX = 200.0

nsteps = int(TMAX / STEP)
x, y = euler_method(T0, C0, STEP, dydx, nsteps)

Hier haben wir die Anfangsbedingungen C0 = np.array([0.0, 0.001, 0.0]) gewählt, was bedeutet, dass zu Beginn nur die Spezies $\mathrm{B}{\mathrm{r}}^{-}$ mit einer Konzentration von $0.001\mathrm{M}$ vorhanden ist. Die Spezies $\mathrm{HBr}{\mathrm{O}}_2$ und $\mathrm{C}{\mathrm{e}}^{4+}$ liegen nicht vor.

Dieser Code-Block schafft es jedoch nicht, uns die richtige Lösung zu liefern. Der Code wird zwar ohne Fehlermeldung ausgeführt, wir erhalten jedoch eine Reihe von Warnungen, wie z.B.

/path/to/your/python_script.py:line: RuntimeWarning: overflow encountered in xxxxx

Sollten Sie das Ergebnis plotten, werden Sie feststellen, dass die Komponenten in y betragsmäßig sehr groß werden, was zu einem arithmetischen Überlauf (engl. overflow) führt. Das kann ein Zeichen dafür sein, dass wir die Schrittweite $h$ zu groß gewählt haben und das Euler-Verfahren instabil wird.

In diesem Fall müssen wir die Schrittweite $h$ tatsächlich auf h = 0.0005 verkleinern, um eine stabile numerische Lösung zu erhalten:

CX_0 = 0.0    # M
CY_0 = 0.001  # M
CZ_0 = 0.0    # M
C0 = np.array([CX_0, CY_0, CZ_0])

T0 = 0.0
STEP = 0.0005
TMAX = 200.0

nsteps = int(TMAX / STEP)
x, y = euler_method(T0, C0, STEP, dydx, nsteps)

Es sei angemerkt, dass sich dadurch natürlich auch die Anzahl der Schritte erhöht; hier waren es $200/0.0005=400.000$ Schritte. Für Anfangsbedingungen $c_Y^0>0.001$ wären die DGLs sogar noch schwieriger zu lösen und wir müssten eine noch kleinere Schrittweite verwenden.

Nun plotten wir das Ergebnis:

c_x, c_y, c_z = y * 1000.0  # convert to mM

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(x, c_x, label='[HBrO2]')
ax.plot(x, c_y, label='[Br-]')
ax.plot(x, c_z, label='2 [Ce4+]')

ax.set_xlabel('time / s')
ax.set_ylabel('concentration / mM')
ax.set_xlim(0, 200)
ax.set_ylim(-0.1, 2.0)

fig.tight_layout()
ax.legend(loc='upper right')

plt.show()

Wir haben hier mit den Funktionen ax.set_xlim und ax.set_ylim die Achsen eingeschränkt. Außerdem haben wir das Argument loc='upper right' an die Funktion ax.legend übergeben, um die Legende nach oben rechts zu verschieben. Das voreingestellte Argument ist loc='best', was die “beste” Position für die Legende automatisch bestimmt.

Wir erhalten das folgende Diagramm: Man erkennt hier den periodischen, impulsartigen Verlauf von $[\mathrm{Z}]$, was in dem obigen Video der Konzentration von Fe(III) entspricht ($[\mathrm{C}{\mathrm{e}}^{4+}]$ in der Originalreaktion) und dort als Blaufärbung zu sehen ist.

Woher weiß man nun, ob die Schrittweite $h$ klein genug ist? Eine Faustregel besagt, dass man bei erfolgter Rechnung mit gegebenem $h$ die Rechnung mit halbierter Schrittweite $h/2$ erneut durchführen soll. Bleibt das Ergebnis gleich wie mit $h$, dann ist $h$ klein genug.

Unabhängig davon, ob $h=0.0005$ klein genug ist oder nicht, konnten wir erkennen, dass das Euler-Verfahren Schwierigkeiten hat das Gleichungssystem (2.8) zu lösen. Wir werden uns im folgenden Abschnitt mit Methoden befassen, die trotz größerer Schrittweiten stabile Lösungen liefern können. Dabei sei daran erinnert, dass das Euler-Verfahren nur den konstanten und linearen Term der Taylor-Entwicklung berücksichtigt (vgl. Gl. (2.5)). Es liegt daher nahe, auch höhere Ordungen zu berücksichitgen, was zu einer Familie von Methoden führt, die als Runge-Kutta-Verfahren bekannt sind.

Übung

Aufgabe 2.1: Lösen des klassischen harmonischen Oszillators mit Euler-Verfahren

Die Bewegungsgleichung eines harmonischen Oszillators ist durch die Differentialgleichung zweiter Ordnung

$$\begin{array}{c}\frac{d^{2}x}{d{t}^2}=-{\omega }^{2}x\end{array}$$

gegeben, wobei $\omega =\sqrt{k/\mu }$. Um diese Gleichung mit dem Euler-Verfahren zu lösen, muss sie zunächst in ein System von Differentialgleichungen erster Ordnung umgeformt werden.

(a) Umformen in ein System von Differentialgleichungen erster Ordnung

Zeigen Sie, dass die obige Differentialgleichung in das System von gekoppelten Differentialgleichungen erster Ordnung

$$\begin{array}{rl}\frac{dx}{dt}& =v\\ \frac{dv}{dt}& =-{\omega }^{2}x\end{array}$$

umgeformt werden kann, indem Sie die Substitution $v=\frac{dx}{dt}$ verwenden.

(b) Implementieren des Euler-Verfahrens

Implementieren Sie das Euler-Verfahren, um das System von Differentialgleichungen erster Ordnung aus Teilaufgabe (a) zu lösen. Gehen Sie dazu wie in der Vorlesung vor, indem Sie die Funktionen dfdt, welche die rechte Seite der Differentialgleichungen berechnet, euler_step, welche einen Schritt des Euler-Verfahrens durchführt, und euler_method, welche das Euler-Verfahren für eine gegebene Anzahl von Schritten durchführt, implementieren.

(c) Herleitung der analytischen Lösung

Zeigen Sie, dass sie analytische Lösung der Bewegungsgleichung des harmonischen Oszillators mit Anfangsbedingungen $x(0)=x_0$ und $v(0)=v_0$ gegeben ist durch

$$\begin{array}{rl}x(t)=& \frac{v_0}{\omega }\sin (\omega t)+x_0\cos (\omega t)\\ v(t)=& v_0\cos (\omega t)-x_0\omega \sin (\omega t).\end{array}$$

Nutzen Sie dazu den allgemeinen Lösungsansatz $x(t)=A\sin (\omega t)+B\cos (\omega t)$.

Plotten Sie die analytische Lösung und die numerische Lösung des Euler-Verfahrens für die Anfangsbedingungen $x_0=1$, $v_0=0$, $\omega =1$ mit einer Schrittweite von $dt=0.1$ für $t\in [0,10]$. Was beobachten Sie?

Runge-Kutta-Verfahren

Die Runge-Kutta-Verfahren sind eine nach Carl Runge und Martin Wilhelm Kutta benannte Familie von Methoden zur numerischen Lösung von Anfangswertproblemen für gewöhnliche Differentialgleichungen. Diese Methoden berechnen iterativ die Lösung im nächsten Zeitschritt aus einer Linearkombination von dem Funktionswert und den Steigungen an verschiedenen Stellen. Sprich man jedoch von dem Runge-Kutta-Verfahren, ist das klassische Runge-Kutta-Verfahren gemeint, was aber nur ein Spezialfall der Runge-Kutta-Verfahren darstellt.

Theoretische Grundlagen

Wir betrachten wieder das Anfangswertproblem (AWP) in Gl. (2.4). Genau so wie beim Euler-Verfahren, wählen wir zunächst ein gleichmäßiges Grid $\{x_n{\}}_{n=0,\cdots ,N}$ mit $x_n=x_0+n\cdot h$ und $h$ als Schrittweite, sowie eine Anfangsbedingung $y_0=y(x_0)$. Dann entwickeln wir $y(x_{n+1})$ wieder in eine Taylor-Reihe um $y(x_n)$: $$y(x_{n+1})=y(x_n)+h{y}^{\prime }(x_n)+\frac{h^2}{2}{y}^{\prime \prime }(x_n)+\frac{h^3}{6}{y}^{\prime \prime \prime }(x_n)+\frac{h^4}{24}{y}^{(4)}(x_n)+\cdots .\text{\hspace{1em}(2.9)}$$ Man spricht von einem Runge-Kutta-Verfahren der (Konsistenz-)Ordnung $p$, wenn das Taylorpolynom bis zum $p$-ten Grad im Verfahren berücksichtigt wird.

Ein allgemeines Runge-Kutta-Verfahren für das AWP in Gl. (2.4) ist dann durch $$y_{n+1}=y_n+h\sum _{i=1}^s{b}_i{k}_i\text{\hspace{1em}(2.10)}$$ gegeben, wobei $$k_i=f(x_n+h{c}_i,y_n+h\sum _{j=1}^s{a}_{ij}{k}_j).$$ Dabei bezeichnet $s$ die Stufe des Verfahrens. Die Koeffizienten $a_{ij}$, $b_i$ und $c_i$ sind charakteristische Parameter des Verfahrens. Hier werden wir nur explizite Runge-Kutta-Verfahren betrachten, bei denen die Koeffizienten $c_1=0$ und $a_{ij}=0$ für $j\ge i$ gelten.

Weil die allgemeine Form doch sehr unhandlich ist, betrachten wir zuerst ein zwei-stufiges Verfahren zweiter Ordnung.

Runge-Kutta-Verfahren zweiter Ordnung (RK2)

Nach Gl. (2.10) können wir ein explizites zwei-stufiges Runge-Kutta-Verfahren zweiter Ordnung formulieren als $$\begin{array}{rl}{y}_{n+1}& =y_n+h(b_1{k}_1+b_2{k}_2)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+h{c}_2,y_n+h{a}_{21}{k}_1).\end{array}\text{\hspace{1em}(2.11)}$$

Für ein Verfahren zweiter Ordnung benötigen wir ein Taylorpolynom zweiten Grades, in dem die zweite Ableitung von $y$ vorkommt. Da wir aber $y$ nicht kennen, haben wir keinen direkten Zugriff auf die zweite Ableitung. Wir können jedoch mit Hilfe der Kettenregel $y^{\prime \prime }$ wie folgt ausdrücken: $$\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{x}^2}=\frac{\mathrm{d}}{\mathrm{d}x}{y}^{\prime }=\frac{\mathrm{d}f(x,y(x))}{\mathrm{d}x}=\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}{y}^{\prime }=\frac{\partial f}{\partial x}+f\frac{\partial f}{\partial y}.$$ Das liefert und das Taylorpolynom zweiten Grades: $$y(x_{n+1})=y(x_n)+hf(x_n,y_n)+\frac{h^2}{2}\left[\frac{\partial f}{\partial x}+f\frac{\partial f}{\partial y}\right](x_n,y_n)+\mathcal{O}(h^3).\text{\hspace{1em}(2.12)}$$

Anschließend entwickeln wir $k_2$ in Gl. (2.11) linear um $(x_n,y_n{)}^{⊺}$: $$k_2=f(x_n,y_n)+h{c}_2\frac{\partial f}{\partial x}(x_n,y_n)+h{a}_{21}{k}_1\frac{\partial f}{\partial y}(x_n,y_n)+\mathcal{O}(h^2).$$ Setzen wir dies in Gl. (2.11) ein, erhalten wir $$\begin{gathered} y_{n+1}=y_n+h(b_1+b_2)f(x_n,y_n)+h^2{b}_2\left[c_2\frac{\partial f}{\partial x}(x_n,y_n)+a_{21}\left(f\frac{\partial f}{\partial y}\right)(x_n,y_n)\right]+\mathcal{O}(h^3). \\ \text{\hspace{1em}(2.13)} \end{gathered}$$

Damit das Verfahren tatsächlich eine Konsistenzordnung von $p=2$ hat, müssen die Koeffizienten vor den Funktionen $f$ und ihren Ableitungen in Gl. (2.12) und (2.13) übereinstimmen, da diese für beliebige Funktionen $f$ gelten müssen. Das führt zu den Bedingungen $$\begin{array}{rl}{b}_1+b_2& =1\\ b_2{c}_2& =\frac{1}{2}\\ b_2{a}_{21}& =\frac{1}{2}.\end{array}\text{\hspace{1em}(2.14)}$$ Das ist ein unterbestimmtes Gleichungssystem mit drei Gleichungen für vier Unbekannte. Wir können demnach einen der Koeffizienten frei wählen und erhalten so eine Familie von konsistenten Runge-Kutta-Verfahren zweiter Ordnung mit zwei Stufen.

Wählen wir $b_1=b_2=\frac{1}{2}$, $c_2=1$ und $a_{21}=1$, erhalten wir das sog. Heun-Verfahren. Wählen wir dagegen $b_1=0$, $b_2=1$, $c_2=\frac{1}{2}$ und $a_{21}=\frac{1}{2}$, erhalten wir das sog. Mittelpunktsverfahren. Neben diesen beiden geläufigen Verfahren kann man natürlich auch andere Kombinationen der Koeffizienten wählen, solange die Bedingungen in Gl. (2.14) erfüllt sind. Eine allgemeine Parametrisierung der Koeffizienten lautet $$\begin{array}{rl}{a}_{21}& =\alpha \\ b_1& =1-\frac{1}{2\alpha }\\ b_2& =\frac{1}{2\alpha }\\ c_2& =\alpha .\end{array}$$

Die Angabe der einzelnen Koeffizienten auf diese Weise ist jedoch nicht sehr übersichtlich, wenn man bedenkt, dass es auch Verfahren mit mehr als zwei Stufen gibt. Daher gibt es eine kompaktere Schreibweise, das sog. Butcher-Tableau.

Butcher-Tableau

Das Butcher-Tableau ist eine übersichtliche Darstellung der Koeffizienten $a_{ij}$, $b_i$ und $c_i$ eines Runge-Kutta-Verfahrens. Das Tableau für ein allgemeines $s$-stufiges Verfahren lautet $$\begin{array}{cc}\vec{c}& \mathbf{A}\\ & {\vec{b}}^{⊺}\end{array}=\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& \cdots & a_{1s}\\ c_2& a_{21}& a_{22}& \cdots & a_{2s}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ c_s& a_{s1}& a_{s2}& \cdots & a_{ss}\\ & b_1& b_2& \cdots & b_s\end{array},\text{\hspace{1em}(2.15)}$$